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ON THE FORMATION OF SHOCK-WAVES IN 
SUPERSONIC GAS JETS 
(TWO-DIMENSIONAL FLOW) 

By D. C. PACK (University College, Dundee) 


[Received 8 September 1947] 


SUMMARY 

When a jet of gas issues from an orifice as a parallel stream with a given super- 
sonic velocity and flows in a steady state through an outer medium at rest, its 
behaviour is governed by the ratio between the exit pressure of the jet and the 
pressure of the outer medium. If this ratio is only a little greater than unity, the 
jet has a periodic structure to a first approximation. This state has been examined 
by earlier workers; it is discussed here from the point of view of the ‘characteristics’ 
of the hyperbolic second-order partial differential equation of potential flow. The 
periodic structure ceases to give an adequate representation of the jet as the pressure 
ratio is increased, and shock-waves occur on account of the compressive effect of 
the outer medium. A method is given for computing the conditions in a steady 
two-dimensional supersonic jet. It is shown how the point of origin of a shock- 
wave and the shape of the shock-wave formation may be obtained by theoretical 
means. The results of two calculations, given graphically, are discussed and com- 
pared, as far as possible, with experimental work. 

An expression is obtained for the minimum pressuré in a jet. 

From the general behaviour of a jet with increasing chamber pressure it is found 
possible to infer the initial direction of the shock-wave arising at or near the muzzle 
of a gun after firing. 


1. Introduction 
THERE are many physical problems which depend for their proper resolu- 
tion upon a knowledge of the behaviour of a jet issuing from an orifice 
with a velocity greater than the local velocity of sound, two of the most 
important arising in connexion with the gas stream from a gun and from 
the exhaust of a rocket motor. These problems are at best rotationally 
symmetrical, and involve, in their theoretical solution, a great amount 
of labour in computation, since analytical solutions proper to the initial 
conditions are not known at present, and numerical (or graphical) methods 
must be employed. The same difficulty arises when two-dimensional 
(plane) jets are considered, but the numerical work is shorter; we shall 
limit ourselves in the present paper to a study of properties of such jets, 
and in particular to those properties which are the most interesting 
features of the natural problems—namely, the formation and the form 
of the shock-waves which occur. 

The existing literature on the theory of supersonic gas jets is concerned 
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with the case in which the pressure of the gas streaming out of the orifice 
is only a little greater than the pressure in the outside medium; the 
solutions of Prandtl (5) represent perturbations about a steady parallel 
jet issuing with pressure equal to that of the outer medium, and the jet 
structure is periodic to a first approximation. As the exit pressure is raised 
the approximations break down, and a more exact examination of the 
gas flow is required. 

On the experimental side, Prandtl (6) took photographs of the plane 
jet issuing from a nozzle designed to give an approximately parallel exit 
stream, for a range of ratios of chamber pressure p, to outside pressure II 


SGREGE 


(ii) (iii) (iv) 
Fie. 1. 


varying from p,/II = 4-16 to p,/II = 8-13; the velocity of efflux, using the 
one-dimensional theory (Reynolds, 7) for the flow in the nozzle of throat 
width 7-6 mm., exit width 11-4 mm., was equal to 1-85 times the local 
velocity of sound. Hartmann and Lazarus (3) carried out experiments 
with a circular nozzle for pressure ratios in the range 3 to 7, the velocity 
of efflux being equal to the local velocity of sound. We note that there 
are three parameters for each jet problem for a given gas and adiabatic 
flow from a reservoir; (1) a scale factor, since jets may be scaled linearly 
with respect to x and y as the width of the orifice is changed; (2) the Mach 
number of efflux; (3) the ratio of exit to outside pressure. The experi- 
mental information at present available does not represent a systematic 
investigation into the effect of the variation of the last two of these 
parameters. 

Although shock-waves are inconspicuous in Prandtl’s photographs, they 
are already strikingly in evidence in those of Hartmann and Lazarus for 
pressure ratios lying well within the range used by Prandtl for his plane 
jets. As p,/II increases (for the cylindrical jet) a short wave occurs in the 
form of a truncated cone without ends, with its rim at some distance from 
the axis of the jet (Fig. 1(i)); the rim moves inwards until a full cone is 
formed (Fig. 1 (ii) ); at still higher pressures a ‘Mach’ or ‘bridge’ wave is 
observed (Fig. 1 (iii) and (iv) ). 

In the present paper the possibility of accounting theoretically for the 
occurrence of shock-waves and of calculating their position has been 
investigated by the solution of two problems in two-dimensional flow using 
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ON THE FORMATION OF SHOCK-WAVES 3 


the method of characteristics. The details of the calculation are given 
in subsequent sections, together with a graphical representation of the 
propagation of the characteristics (‘wavelets’) through the flow, and the 
positions of the shock-waves. The waves of rarefaction which set in at 
the bounding edges of the orifice are reflected from the jet boundary as 
waves of compression, and the diagrams show how this compression 
involves a sharply steepening pressure gradient; the shock-wave begins 
where the gradient becomes infinite. 

The exit velocity of the jet in each case corresponded to a Mach number 
of 1-5, and the ratios of chamber pressure to external pressure were 7-1 
and 20 respectively. The former, lying within the range of Prandtl’s 
experiments, gave rise to a shock-wave of type corresponding to that in 
Fig. 1 (i); the latter resulted in a shock-wave which began as in Fig. 1 (iv), 
but which crossed the axis instead of giving a Mach wave. There seems 
no reason to suppose that such a wave might not be observed experi- 
mentally under appropriate conditions; the pressure ratio was well above 
those used in the experiments quoted. 


2. The equations governing the flow 

In order to examine the properties of the gas motion, we consider the 
jet to have reached a steady state, and to be moving irrotationally, so 
that a velocity potential ¢ exists. This means that we neglect any transfer 
of vorticity which may occur as the result of the discontinuity of velocity 
which exists at the boundary between the jet and the outer medium. 

Let axes of coordinates x, y be taken in a normal cross-section of the 
jet, the x-axis lying along the axis of the jet, 
and the y-axis perpendicularly to it. Let the 


origin be taken in the orifice (Fig. 2). Let pp, — a 





Po; & be the pressure, density, and velocity of HH, o4— 

° . ° —— x 
sound, respectively, of the gas at rest in the — 
discharging vessel (reservoir); write a)u, av 


, ; : Fic. 2. 
for the velocity components in the (x, y)-direc- 


tions respectively, and a,a for the local velocity of sound, so that wu, v, a 
are non-dimensional quantities. If q? = u?+v?, so that q is the non- 
dimensional resultant velocity, and if @ denotes the angle between the 
direction of motion at any point and the z-axis, then tan@ = v/u, and we 
may define the Mach angle m = sin-1(a/q). 

Let the gas obey the adiabatic law so that the density p and the pressure 
p are connected by the equation pp-” = constant, where y is the ratio 
of the specific heats. If we write P = p/p), making P a non-dimensional 
pressure symbol, then P = a?vy-, 
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The equation for the velocity potential ¢ is 


a9 
eo 
(u?—a?) ih + 2uv 


ox" OX 


a2 
9\ OC 

FE + (ota ¥ = 0 (1) 
oy dy? 

for two-dimensional flow; since the stream-lines issue from a region of 
constant conditions (the discharging vessel), Bernoulli’s equation possesses 


a universal constant for the flow, and gives the relation between the gas 


a 


velocity q and the sound velocity a: 

gq? = 2(1—a?)/(y—1). (2) 
If the half-width of the jet be taken as the unit length, it is clear from 
(1) that the results may be scaled by multiplying x and y by any one and 
the same factor. 

The gas flow may be followed by considering the changes which occur 
along the ‘characteristics’ (see, for example, Goursat, 2) of the differential 
equation (1). The differential equation of the characteristics is 

(u?—a?)u?— 2uvp+ (v?—a?) = 0, (3) 
where p = dy/dz. 

At supersonic velocities the characteristics are real, giving two families 


of lines; one characteristic of each family passes through each point of 


the flow. 
If the roots of (3) are uw, and pg, it is found that: 


—— , pe ee ae 
along the family of characteristic lines Fe Hy 
Ax 
we have du+p.dv = 0 
eee ‘ ad eas dy 
and along the family of characteristic lines . = ple; . (4) 
dx 


we have du+p,dv = 0 
It is easy to show that the stream-line through a point bisects the angle 
between the characteristics through the point, and that the characteristics 
make the Mach angle m with the direction of flow. 
Using the definition of m and the equation (2), we may write g = q(m), and 


hence we have , 
. u = q(m)cos 8, v = q(m)sin 8, 


fy = tan(é+m), fl, = tan(@—m). 
Changing the variables in the relations (4) to 6, m, we obtain integral 
relations along the characteristics: 


along 3 = tan(@+m), p(m)+6 = constant, (5) 
dx 

dy : 

and along - = tan(@—m), p(m)—é = constant, (6) 
dx 


where p(m) = xtan(; tan m)—m, and A? = (y—1)/(y+1). 
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It should, perhaps, be pointed out that the expressions on the left-hand 
sides of these conditions are constant, in general, only along individual 
members of the family. An exception occurs when a set of characteristics 
issues from a region of uniform flow (or of rest). Relations similar to (5) and 
(6) were first found by Steichen (8) in his dissertation on two-dimensional 
gas motion. 


3. The properties of a plane supersonic gas jet 

The behaviour of the gas jet issuing from an orifice at supersonic velocity 
into a medium of lower pressure, in the absence of shock-waves, may be 
followed from Fig. 3, bearing in mind the relations along the charac- 
teristics. Let M, M’ be the corners of the (section of the) orifice, p, the 


o-TT 





pressure across MM’, and II the external pressure. The stream flows 
unchanged in the region A, M’, in which the characteristics form two 
sets of parallel lines making the Mach angle with the axis of flow. We 
shall refer to these characteristics as being of a-type; coming from a region 
of uniform axial flow, they possess the same constant in (5) and (6). 
Around M and M’ the gas expands freely in the regions MA, B, M’A, B’, 
and then flows in parallel lines in the regions M BC, M’ B’C’ at pressure II. 
Characteristics spring from the corners in the two regions of expansion; 
these we call characteristics of B-type, each possessing its own particular 
‘constant’, defined by equation (5) or (6) together with the well-known 
expansion round a corner (‘Meyer expansion’—see, for example, Taylor 
and Maccoll, 9). Since there is a universal constant for either set of 
x-characteristics, it follows that when such a set is met by a f-charac- 
teristic, there is a unique @ and a unique m holding. Consequently, under 
such circumstances, a 8-characteristic is straight, and the a-characteristics 
intersect it at a constant angle. This is the case in the regions of the 
Meyer expansion MA, B, M’A, B’, and in the regions BA, DC, B’A, D'C’. 
It should be noted that each family of characteristics contains charac- 
teristics of both a- and f-types. 
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Across A,D, A, D’ a-characteristics pass, making equal and opposite 
angles with the axis of the jet. It is easily seen that this implies a rhom- 
boidal region A, DA, D’ in which the pressure is uniform (equal to ps, say), 
and in which the gas streams parallel to the axis. 

The jet boundary is a sheet of velocity discontinuity, with pressure 
constant across it, and it reflects a characteristic into the other family 
(see Fig. 3), the angle of reflection being equal to the angle of incidence 
since the boundary is a stream-line. 

The flow may be solved analytically except where B-type characteristics 
intersect, e.g. A, BA,, CDE. In such regions numerical methods must 
be employed. This also entails numerical methods for the determination 
of points on (for example) CD, which are obtained from BA, via the 
straight-line £-characteristics in BA, DC, using conditions of similarity. 

From the plane DD’ the earlier part of the phenomenon is repeated 
but with the stages in the reverse order, the expansion regions becoming 
regions of compression. The f-characteristics reflected from CH, C’E’ are 
straight in the compression regions DE F'A,, D’ E' F'A, and the compression 
waves interfere in the region A, FA, F’. The conditions in EFG, EF’ F’G’ 
are uniform at pressure II, HG and E’G’ being the stream-lines which form 
the boundary of the jet. 

In the Prandtl-Busemann solution of the jet problem (see Busemann, 1), 
the f-characteristics, after traversing the interference region A, FA, F’, 
become straight on meeting the a-characteristics reflected from EG or 
E’G'’, and meet in a point at G (for one set) and G’ (for the other set). 
The whole phenomenon is then repeated with periodic structure, the 
wave-length being OA;. This implies the reversibility of the flow; thus, 
the flow must be isentropic, and no shock-waves occur. In practice, how- 
ever, as we shall show, f-characteristics of the same family intersect, and, 
as has been indicated by Taylor and Maccoll (9), this may be interpreted 
to imply the existence of a shock-wave. For ratios of exit to outer pressure 
close to unity the Prandtl-Busemann solution is, nevertheless, an ex- 
tremely good approximation, because, in this case, the shock-waves are 
formed in the neighbourhood of G and G@’, and are not strong; and it may 
be shown that, for weak shock-waves, the change in the value of p/p” is 
of the third order in the density change—that is to say, the flow is 
reversible to the third order of small quantities. 

The creation of a curved shock-wave in the interior of a fluid involves 
the breakdown of the irrotational nature of the flow. The increase in 
entropy in crossing the wave is accompanied by the formation of vorticity 
behind it. This appears at first sight to contradict the principle of con- 
servation of irrotational flow, but the latter holds only so long as the flow 
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remains continuous, and it is well known that a compressive motion cannot 
usually be maintained for more than a finite length of time without the 
occurrence of a discontinuity. 


4. The choice of theoretical problems 

In order to see to what extent it might be possible to account theoretically 
for the occurrence and the form of the shock-waves in a supersonic gas jet, 
it was necessary to choose for consideration problems which would not 
involve a prohibitive amount of calculation, but which might be expected 
to provide interesting illustrations. The exit velocity of the parallel stream 





when the 
exit velocity is sonic or near-sonic, difficulties and inaccuracies arise near 
the orifice on account of the steep gradients of the characteristics—the 
ratio of specific heats was taken to be 1-3 (corresponding approximately 
to the gases issuing from the muzzle of a gun), and it was left to decide 
the chamber pressure. For the first problem the ratio of exit pressure to 


at the orifice was taken to correspond to a Mach number of 1-5 


external pressure (p,/Il) was taken to be 2, corresponding to a chamber 
pressure of 7-111; this served to establish the computational procedure. 
For the second problem a chamber pressure of 20II was chosen (p, = 7-5II), 
in order to seek a different type of shock-wave formation. 


5. The method 1 computation 

Using the non-dimensional quantities defined earlier, and a suffix z to 
represent the state in the uniform regions MBC, M’B'C’, the conditions 
required to begin the calculations are as set out in the following table: 





“Pe eee : “cee 
Problem | ay, | Uy Qi |Af| mf P, | ay | In | A m?, Fr. 
lx 














I | 0°864675 | 1°297013 | oO | 41°810/} 0°283612 | 0°798213 555326 | 14°382 30°878 | 0-141806 








II om ” ” ” - 0°707752 | 1°824074 | 32°627 | 22°830 | 0*050000 





It is seen at once, by evaluating the differences of (@—m), that the 
angles of expansion A, MB, A, M'B’ are equal to 25-314 degrees and 
51-607 degrees—and by evaluating the differences of (9+m) it is seen 
that the characteristics M'A, B, MA, B’ are deflected through angles of 
3-450 degrees and 13-647 degrees in traversing the region of expansion— 
in the two cases respectively. 

The characteristic M’A,B in the expansion region A, B follows the 
curve 7?(sin 8)"“"cos 8 = const., where 


B tan-{ tan m) and A? = (y—1)/(y+1) 


as before. The equation may be derived from first principles, using 
the conditions of the Meyer expansion; r represents the distance from 
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the appropriate corner. The constant of the a-characteristics is to be 
determined, and using this in conjunction with the equation of the charac- 
teristic found above, values of a, y, 6, and m are calculated at points on 
A, B corresponding to chosen intervals of the angle of 
expansion, the f-characteristics in this region being 
straight lines through the corner. We make use of the 
axis of symmetry, so that only one-half of the field 
has to be investigated; we restrict ourselves to the upper 
half. Given the values defining the state at a number 
of points on A, B, conditions in A, BA, are calculated, 
using the method of characteristics. The field has to be 
built up step by step, the conditions at two points 7 
and B leading to the knowledge of the state at a new point P which is 
the intersection of the characteristics shown in Fig. 4. 

From (5) and (6), using suffixes to denote the points to which particular 
quantities are referred, we have for the point P: 





Op = £(87+0p)—3(Pr—Pp), (7) 


Pp = 3(Prt+Ppz)—3(I7—4Op). (8) 

The value of mp is deduced from (8) by the aid of previously prepared 

tables of p(m) as a function of m. When the point P lies on the axis of z, 
the formulae reduce to 

Op = 0, (7a) 


Pp = Pr—Or. (8a) 


These formulae are exact. It is necessary, however, to introduce an 
approximation in order to obtain the coordinates of P. The differential 
equations for the characteristics are replaced by finite difference relations, 
and we use, with a slight change of notation, the formulae previously 


established by Hill and Pack (4). If we write u for the gradient of a 


characteristic of the family » = tan(@+m), and mi for the gradient of a 
characteristic of the family » = tan(@—m), always attaching the appro- 


priate suffix, then putting wpt+pp, = ©, upt+pr = &, and replacing all 
gradients by the mean of the gradients at the ends of the are under 
consideration, the coordinates x, y of P may be expressed by the equa- 
tions 


(2—2X)xp = LA _pt(—L)ap+2(yp—yz) (9) 


(X—Z)yp = Ly7+(—Z)ygt+(—4$z X)(xp—zg) 
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When P is an axial point we have 
9 
“YT 


Lp Lao 
(—2) 


(9a) 


Y p 0 


Using the above equations, the region A, BA, is completed, and it 
becomes necessary to calculate points on CD corresponding to those found 
on BA,, remembering that the f-characteristics are straight in BCDA,, 
and that the a-characteristics cross them at constant angles along their 
lengths. 

We may calculate the coordinates (X»,Y,) of C, since we know the 
equations of MC and BC; with these and the 


coordinates of two successive points B(x», y,) and Cn 
B,(x1,y,) on BA,, we may find the coordinates Cres 
(X,,Y,) of the point C, where the f-characteristic 
through B, meets CD. In general (Fig. 5), know- 
ing the coordinates of C,,(X,,,Y,,) on CD, and the p 
consecutive points B,(x,,y,) and B, 44(@n415Yn41) 3 Bn+i 
on BA,, we may calculate the coordinates of Fic. 5. 
C.44(X 44; ¥n4,) by means of the formulae: 
Ynsa—Yn = (Xnsr—XnYnsa—Yn)/ (Cnn) (10) 
Fn+i—Ynsa (Xn sins )Pn+y (11) 


where ,,,, is the gradient of the upward f-characteristic through B,,,,. 
(11) is the equation of this characteristic, and (10) is a line C,,C,,,, parallel 
to B, B,,,, the true gradients of these lines being replaced by the mean 
value of the gradients at the extremities. 

The solution of (10) and (11) gives 


B “ ae (Y,, Pn+1Xn)—(Ynii—PntiTn+1) (12) 


nT ? 


Pnsit4(—2) 


X)i(Yn — Pn 1Xn)—(Yn +1 Fn +17 n a} (13) 
nat $(—2) 


Having obtained conditions at points along CD, we proceed to compute 
values of 6, m for points in CDE, taking into account the fact that the 
pressure is constant on the boundary. Constant pressure implies constant 
m. The coordinates of the points on the boundary may be calculated 
using a slight modification of (9). 


o- 
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Let 6@,,, 6,4, be the directions of the stream at two successive points 
which we seek on the boundary (Fig. 6). Let the coordinates of the nth 
boundary point be used as those of 7’ in the earlier work. Put 

o = tan@,+tan0,,, 


(0,,., 18 known, since m is given and we know the @, m relation on the 
characteristic BP). 





The formulae for the coordinates of the (n+-1)th boundary point P are 


then identical with those of (9), provided we replace X everywhere by oc. 
The boundary curves after the manner of CE in Fig. 3. 

On completion of this part of the work, the state of the flow is known 
at a number of points on DE, and these serve to give the field in the 
region DEF A,, except in so far as this region may be affected by the 
occurrence of a shock-wave; failing the latter, equations similar to (12) 
and (13) may be set up to find the flow conditions on A, F, and the 
computation continues. ! 


6. The origin of a shock-wave 

We have already seen that the effect of the interaction of the expansion 
waves from the opposite faces of the orifice is to increase the expansion 
ofeach. In consequence, the compressions which result from their reflection 
at the free boundary of the jet are so much the greater, and indeed, the 
convergence of the B-type characteristics which represent the compression 
is such that consecutive characteristics eventually intersect one another, 
and the flow field ceases to be single-valued in the quantities expressing 
the state at a given point. For p, ~ 7II, this intersection of characteristics 
of the same family begins in the region A, F'A,, after the opposing com- 
pressions have experienced some interaction; for p, = 20II, the intersection 
occurs between C and D, i.e. before the reflected characteristics reach the 
axis. 


The intersection of characteristics of the same family is to be interpreted 
as the sign of the breakdown of continuous motion and the beginning of 
a shock-wave. Taylor and Maccoll (9), investigating the flow past a thin 
curved surface placed edge-on to a supersonic stream, imagine the shock- 
wave formed as being built up by the intersections of characteristics of the 
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same family, the shape of the shock-wave, to a first approxiniation, being 
the envelope of these characteristics. In fact, the characteristics will only 
be tangential to the shock-wave at its point of origin, where it is of 
infinitesimal amplitude. As the strength of the shock-wave increases, the 
characteristics meet the wave at an angle (see Figs. 7 and 8) determined, 
along with the position of the shock-wave, by an application of the 
equations of flow through an oblique discontinuity. This method will be 
available so long as the flow behind the shock-wave remains supersonic ; 
it will break down at subsonic velocities, because the characteristics then 


become imaginary. 


7. Determination of points on the shock-wave 
The shock-wave begins as a wave of infinitesimal amplitude at a point 
where two characteristics (including CD or A,F’), whose angles of 


FIG. 


departure from M’ differ by an infinitesimal quantity, run together tan- 
gentially. In order to find this point the field is computed without respect 
to intersections. The angles of intersection of the ‘critical line’ (CD or 
A, F) and the successive f-characteristics are calculated, and the point 
at which a characteristic intersects the critical line at zero angle is then 
computed by the method of divided differences. The initial direction of the 
shock-wave is the direction of the critical line at this point. The calculated 
flow downstream of this position must now be abandoned, and the field 
re-evaluated allowing for the presence of the shock-wave. The results of 
the two calculations are given graphically in Figs. 7 and 8, showing the 
shock-waves and some of the characteristics. As the shock-wave curves, 
the motion behind it ceases to be irrotational, but as the curvature is not 
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very great, and our primary interest is qualitative, the vorticity has been 
neglected in the analysis and computation. 

The position of any point P on the shock-wave will be determined (i) by 
the conditions on the characteristics which carry the flow into it; these, 
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being in each case straight lines, can be written in the form y = bate, 
and can be represented parametrically by their appropriate value of m; 
(ii) by the conditions holding on the characteristics which run into what 
may be termed the rear of the wave (the convex side). The problem thus 
reduces to the following (Fig. 9): 

Given a point P of the shock-wave, we have to determine an are PP’ 
of it such that the gas flowing into P’ is deflected in a manner consistent 
with the equations for the passage of gas through an oblique shock-wave, 
and also with the conditions on that characteristic in the rear of the wave 
which meets the latter at P’. 

Let suffixes F, R refer to the front and rear of the shock-wave respec- 
tively. Let %, x be the angles between the deflected and the incident flow, 
respectively, and the normal to the shock-wave at P (Fig. 10). 

The relations connecting the conditions on the two sides of the shock- 
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wave follow from the well-known-equations (see, for example, Taylor and 
Maccoll, 9). If z = p,/p,, where p, is the pressure in front, and p, the 
pressure in the rear of the shock-wave, then 
cos?y {(y—1)+(y+ 1)z}sin’m,, 
tan %/tan x {(y—1) + (y+ L2y/{(y+1I)+(y—1)z}, 
and SIN Mp = sin m,.f(z)cos p/cos x, 


where f(z) = [e{(y—V) + (y+ Dah /{(y+ I) + (y—1)2}}!- 
By means of the above equations it is possible to construct universal 
tables for the greater part of the problem; one needs to tabulate against 
my: (i) x as a function of (’—y); (ii) mp as a function of (fJ— x). Tables 
b, ¢ against m, must be prepared for each problem from the 
characteristic relationships in front of the wave. 

The problem of finding P’ is solved, in practice, by moving along a rear 
characteristic, guessing the value of m, for the point at which it will meet 
the shock-wave, thus fixing @,,, and finding a deflexion of the stream giving 
values of @,, and m, which are consistent with the relation holding on the 
rear characteristic. For this deflexion we then have a determinable value 
of the gradient of the shock-wave, and we must next, by an approximation 
to the equation of the are of the shock-wave, find its intersection with 
the forward characteristic (m,). This gives values for the coordinates 
of P’; it must be checked if these values are consistent with the (approxi- 
mate) equation of the rear characteristic—if not, a new m, must be 
guessed and the procedure repeated. The method simplifies when the 
field in front of the shock-wave is uniform, since m, is fixed. One 
calculation is sufficient for each point. The method would have increased 
difficulty in a rotationally symmetrical jet problem, since the forward 
characteristics would no longer be linear; and indeed, difficulties would 
arise for the same reason, if we were to try to continue the second of the 
problems of this investigation beyond the reflection of the shock-wave at 


the axis. 


8. Discussion of results 

The accuracy of the solutions obtained is difficult to assess, since the 
errors involved in each step are additive; these errors, however, arise only 
in the location of the coordinates. In the present calculations, on the 
basis of the sizes of the blocks of points requiring numerical evaluation, 
and of the differences found by doubling the interval between successive 
points of the network, it seems likely that, beginning with four places of 
decimals in x, y, the positions are correct to two places of decimals, the 
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reliability falling off in those parts of the field which are most remote from 
the orifice, i.e. behind the shock-wave. The results of the two calculations 
are given graphically in Figs. 7 and 8. 

For the jet of Fig. 7 the distribution of pressure along lines parallel to 
the axis of the jet is given in Fig. 11. It may be seen at once how the 
pressure falls rapidly in the expansion region to a value p, which differs 
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but little from Prandtl’s approximate value (II/p,)I], derived in the next 
section. In the region of ‘interference of the two expansion waves there 
is a sudden increase in the rate of fall of pressure, shown clearly by the 
discontinuities in gradient on the graphs for y = 0-25, y = 0-50, at S,, 7, 
respectively. All the lines along which pressure has been calculated pass 
through the region of converging straight characteristics, except for the 
axis which goes directly from the uniform region of pressure p, into the 
part of the field in which the waves of compression reinforce one another. 
This reinforcement of the compression shows itself for all the lines by a 
sudden increase in gradient of the pressure graph (R3,S,, for example), 
and the series of graphs'drawn illustrates clearly the way in which the 
shock-wave begins. The gradient discontinuity originating on the axis at 
R, in Fig. 11 (ie. A, in Fig. 3) becomes more and more severe forwards 
and away from the axis, until eventually it becomes infinite with a result- | 
ing instantaneous rise in the magnitude of the pressure itself, ie. with | 
the formation of a shock-wave. This is almost the case for y = 0-28 in 
the figure. For y = 0-5 there is a finite jump in pressure (at 7, in Fig. 11), 
since this line actually passes through the shock-wave itself. 

For the 2011 chamber-pressure problem (Fig. 8), the shock-wave curves 
concavely towards the axis, where it is reflected. The Mach numbers of 
the flow in front, between, and behind the shock-waves, on the axis, are 
4:17, 3-25, and 2-55 respectively. It is to be observed, in this solution, 
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how sharply the boundary of the jet curves inwards. The jet has almost 
reached its narrowest section at a point where, with periodic structure, 
it would be at its maximum width. 

From the point of view of comparison with experimental work, it is 
to be noted that the lesser chamber pressure gives rise to a shock-wave 
formation of the type shown in Fig. 1 (i), and this accords with expectation, 
taking into account the difference between the experimental results for 
two-dimensional and rotationally symmetrical jets. At the higher pressure 
ratio, the theoretical shock-wave pattern does not conform to the ‘Mach’ 
or ‘bridge’ type, but crosses the axis of symmetry in a regular manner. 
This may be contrasted with the results of the experiments upon a 
rotationally symmetrical jet with sonic efflux velocity, for which the Mach 
wave pattern occurred for all pressure ratios p,/II > 4-3. 

It is to be noted, finally, that with the increase of chamber pressure, 
the point of onset of the shock-wave moves backwards along the critical 
line CD (or C’D’) towards C (or C’). 

9. The minimum pressure 

The value of the pressure p, in the region A, DA, D’ (the minimum 
pressure in the first period) may be deduced from the characteristic 
relationships. The condition on M’A, BC, namely, p(m)+6 = p(m,), 
enables the condition on MBA, to be written: 

p(m)—@ = 2p(m,)—p(m,), 

where the suffix unity refers to conditions on efflux, and the suffix 7 to 
the uniform regions MBC, M’ B’'C’ in which the pressure is equal to that 
of the outer medium. If the suffix be used for quantities in the region 
A, DA, D’, then putting 6 = 0 above, 

p(m,) = 2p(m,)—p(m,). (14) 
Now let an angle v be defined by the equation tanm = Atanv. Let 
p, = nil; then combining the pressure relationship p/p, = a?”y-) with 
Bernoulli’s equation, one obtains 

p/Po {2/(y+1)}7y-V(sin v)2vy-D, 


leading to the two equations which together with (14) determine p,: 


sinv, = (1/n)’-ysin v,, (15) 
sin vg\2vky-1) sin v,\?7ky-)) 
Pe nil| a =n" ; (16) 
a sin v, \sin v,, 


SiN Vp 


. s 2yKy—l) 7 . 
We see that p,p. = n2TT°( . If we put = l-+e, where e« is 


sinv, 
small compared with unity, it is easy to show that (sin v,/sin v,)?”7-) ~ n-2, 
We thus derive the well-known Prandtl rule for the minimum pressure, 
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P,P, ~ Il?. The attainment of the minimum pressure expected from the 
equations (14)-(16) is clearly dependent upon the flow in MOM’A, being 
totally unaffected by the presence of any shock-wave; shock-waves-may 
freely exist outside this region without affecting the result. 

It has already been remarked that the minimum pressure found with 
Pp, = 211 was in fairly close agreement with that given by Prandtl’s rule. 
For the pressure p, = 7-5I1, however, the minimum pressure was calculated 
to be 0-077II, differing greatly, as we should expect, from the value derived 
from the simple rule. 


10. The flow of gases from the muzzle of a gun 

It was noted at the end of §8 that the effect of an increase in the 
chamber pressure is to move the point of onset of the shock-wave back 
along G’C towards C. The point cannot move beyond C because it must 
come in a region where the motion is one of compression; on the other 
hand, C can be made to approach M by reducing the velocity of flow 
of the gas from the orifice. Using these facts as a starting-point we can 
infer that, for sonic exit velocity and a sufficiently high chamber pressure, 
the shock-wave may be considered to begin, with zero strength, approxi- 
mately at the orifice itself. 

For a parallel, uniform jet issuing with a pressure in excess of that of 
the external medium, this seems to be the only way in which a shock- 
wave can begin at the orifice. For the same jet issuing at under-pressure 
a shock-wave beginning at the lip of the orifice is to be expected; and it 
is possible for shock-waves to occur inside the lip of a nozzle for which 
the design is not sufficiently accurate to ensure a parallel uniform stream. 

As a corollary to the above remarks, it is of interest to consider the 
gas flow issuing from the muzzle of a gun; the flow velocity is approxi- 
mately sonic. Since in the limit at M the initial directions of the charac- 
teristics in the two-dimensional and the rotationally symmetric problems 
are the same, we may infer, in particular, that the direction of the 
characteristic CD at C in the plane case, in the limit as sinm,— 1 and 
C approaches M, gives the initial direction, approximately, of the shock- 
wave from the muzzle of a gun. 

Using the notation of the earlier part of the paper, it follows from 
Bernoulli’s equation that 

{3 cosec’m,,+- 1/(y—1) }19-Y/v = {4 cosec?m,+-1/(y—1)}(nI1)v-Y"”, 


> 


BC makes an angle (6,+-m_) with Ox, and it is known that 


—6,, = p(m,)—p(m,). 


The angle between MC and Ox is 6,, hence the direction at C of CD, the 
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reflection of BC, is at an angle (m,—6,) with Oz, i.e. it makes an angle 
& = 4n—m,+ 6, with the direction of MM’, which, by the relation above, 


gives ¢ ka—p(m,)—m,,+p(m,). For sonic efflux, m, = $a and 
p(m,) }ar(1—A)/A. 
From Bernoulli’s equation, tanm, = A/(nv-Y/y—1)*. It follows at once 
that 
3 x tan ty (n'y—Mly— 1). (17) 
For n = ©, this gives € = 7/2A ~ 249°; the shock-wave begins in a 


direction inclined steeply backwards from the plane of the orifice. 
For n = 500, a high muzzle pressure for a gun,  ~ 168°, i.e. it begins 
just forward of the continuation of the plane of the muzzle. 
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A VARIATIONAL PRINCIPLE OF MAXIMUM 
PLASTIC WORK IN CLASSICAL PLASTICITY 


By R. HILL (Cavendish Laboratory, Cambridge) 
[Received 8 October 1947] 


SUMMARY 

The classical equations of Lévy-Mises and Prandtl-Reuss for an ideally plastic 
material are reviewed. A variational principle of maximum plastic work is derived 
for plastic states of stress satisfying the Lévy-Mises relation and the Huber-Mises 
yield criterion. Uniqueness theorems are established for a completely plastic body 
under prescribed boundary conditions. The variational principle is applied to the 
problem of a uniform bar of arbitrary section deformed in combined tension, torsion, 
and bending. 


1. The yield criterion and stress-strain relations 

THE following considerations are concerned with an ideally plastic material 
which is perfectly elastic up to a sharp yield point and thereafter does 
not work-harden. The criterion of yielding under combined stresses o;; is 
taken to be that of Huber-Mises: 


05; 0%; = 2k3, (1) 


where o;; is the deviatoric component of the stress tensor a;;. o;; is also 
called the reduced stress tensor. k is equal to Y/V3, where Y is the yield 
stress in uniaxial tension or compression. Since the work-hardening is 
assumed to be zero, k remains constant during the plastic deformation. 
The relations between stress and strain-increment for an isotropic 
element of material which is being plastically deformed are taken to be 





1— 2p 
and de,;; = Ev (3) 


E, G, v are the elastic constants: Young’s modulus, shear modulus, and 
Poisson’s ratio respectively. de;; is the tensor representing an increment 
of true or natural strain, measured with respect to the current configura- 
tion. dA is a scalar, essentially positive, but otherwise unspecified and 
generally varying both in time and space. The first term on the right-hand 
side of equation (2) is the deviatoric component of the increment of elastic 


strain, also measured with respect to the current configuration. The 
second term is the increment of plastic or permanent strain. Equation (3) 
expresses the experimental fact that the elastic compressibility is un- 
changed by the plastic distortion. It is supposed, further, that all the 
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elastic constants remain invariable, provided they are defined in terms 
of the current configuration. These equations are due to Reuss (13), who 
based them on the work of Saint-Venant, Lévy, von Mises, and Prandtl. 
When the plastic flow is ‘free’ and takes place under constant stress, the 
elastic strain-increments are zero and the equations are then equivalent 
to the Lévy-Mises relations: 
de;; = o;;dA. (4) 

The criterion of yielding (1) has been shown to be valid for copper, 
nickel, aluminium, iron, strain-hardened mild steel, medium carbon and 
alloy steels (Ros and Eichinger, 14; Taylor and Quinney, 19; Davis, 1 
and 2; Lessells and MacGregor, 6). Of course for a real metal k and Y 
depend on the previous strain-history. The Lévy-Mises relations (4) have 
been shown to be valid to a first approximation for copper, nickel, iron, 
aluminium, and strain-hardened mild steel (Lode, 8; Taylor and Quinney, 
19: Schmidt, 17; Davis, 2). When the flow is not free, and the elastic 
and plastic components of the strain are comparable, there is enly limited 
evidence in support of the Reuss equations (2). The general appropriate- 
ness of the Prandtl-Reuss method of allowing for varying elastic strains 
in the ideally plastic body can, however, hardly be doubted. 
2. Stationary work principle and the plastic potential 

A method of representing the reduced stress and strain-increment 
tensors on a plane diagram was introduced by Reuss (13). This is parti- 
cularly useful for a general discussion of stress-strain relationships. The 
principal components of stress o,, o), o3 are plotted along Cartesian axes 
of reference (Fig. 1). OS is the vector (0,,0¢,,0,) and OP the vector 
representing the reduced principal stresses 
(o;,0,03). P therefore lies in the plane S(0;,02.05) 


Rove. ) 


Fic. 1. Reuss’s geometrical repre- 
sentation of stress. 


o,+0,+0, = 0. The vector PS is perpen- 





dicular to this plane, and has components 
(o,0,0), where o 4(o,-+0,-+03) is the 
mean of the principal stresses. The prin- 
cipal components de,, de,, de, of the incre- 
ment of true strain can be represented bya _ ig the plane o,-+0,+-0, = 0. 
vector in the same space by plotting 2G de,, 

ete. The vector representing the reduced strain-increment (de, deg, des) 
will lie in the plane o,+0,-++-0, = 0. In this way the relations between 
reduced stress and strain-increment for an isotropic material can be 
followed on a plane diagram; the hydrostatic component of stress need 
not be considered in discussing plastic behaviour. The strain-history of 
an element can be represented by some curved path in the plane, built 
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up of successive vector increments. The ‘final’ strain will be denoted by 


the point Q. There is no physical significance in the vector 0Q, except 
for purely elastic strains or when the strain path is such that the principal 
axes of strain-increment do not rotate relatively to the element. 

Suppose now that the principal axes of stress do not rotate relatively to 
the element. While the element is still only strained elastically, P and 
Q coincide. ‘After yielding, the points separate and P moves on a circle, 
centre O and radius ,/(3)Y, representing the yield criterion (1) (Fig. 2). 


a} 
nef 
Q 
/ 





Fic. 2. Huber-Mises circle with stress and strain-increment 
vectors in plane o,+-0,-+0; = 9. 
Suppose QQ’ is a strain-increment given to the element following a 
moment when the stress vector is OP, and suppose the stress changes to a 
neighbouring stress vector OP’ on the circle during this strain. The principal 
axes of the elastic and plastic strain-increments coincide in view of the 
initial supposition. According to equation (2) the plastic strain-increment is 
parallel to OP, and the elastic strain-increment is parallel to PP’, that is, 
perpendicular to OP. Hence, if RQ’, QR are the components of QQ’ along 
and perpendicular to OP, they represent respectively the plastic and elastic 
strain-increments. The change in stress PP’ is thus equal to QR. The 
work done on the element is positive so long as OP.QQ’ is positive; dA is 





then also positive. Otherwise unloading occurs. 

Consider, now, free flow where the elastic strain-increment is zero or 
negligible. We do not need to retain the restriction that the principal 
stress axes do not rotate. Suppose the strain-increment QQ’ is prescribed, 
while no restriction is placed on the stress beyond requiring P to lie on 
the circle. Then the incremental work done in the strain QQ’, correspond- 


, 7 QQ’ —_ —~ —— ms: 
ing to some stress P, is just ee the projection of OP along QQ’. This 


is a maximum when OP is parallel to QQ’, that is, when the Lévy-Mises 
relation is satisfied. Thus the Lévy-Mises relation corresponds to maxi- 
mum plastic work in a prescribed strain-increment for varying stress 
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subject to the Huber-Mises yield criterion (Hill, Lee, and Tupper, 4). 
Von Mises (10) stated this merely as a stationary work principle, without 
examining whether a maximum or minimum was involved. He suggested 
further that for any yield criterion the appropriate stress-strain relations 
could be generated on the basis of such a principle; this idea has also been 
independently proposed by Taylor (18). The stress corresponding to 
a given strain-increment would then be represented by the point on the 
yield locus where the normal is parallel to the increment of strain vector. 
In 1928 von Mises introduced the concept of a plastic potential. If the 
criterion of yielding is 
f(2;;) - constant, (5) 
the plastic potential is defined to be the function f of the stress com- 
ponents. The definition is not quite unique, but this does not matter in 
the subsequent application. Von Mises suggested that the stress-strain 
relations in free flow which are appropriate for a material with this yield 
criterion should be 
de,; = Ft a, (6) 
C0;; 
where dA is a positive scalar factor of proportionality. Since o;,; is a tensor, 
so also is éf/éo;;; and if f is invariant in form for an arbitrary choice of 
axes (as it must be for an isotropic material) then so is @f/é0;;.. Hence the 
above relation is thus far suitable to represent a strain. When f is the 
Huber-Mises expression o;;0;; it is easily verified that equation (6) leads 
to the Lévy-Mises relations. In general, as von Mises showed, postulating 
the existence of a plastic potential is equivalent to postulating a stationary 


work principle. For the condition that o,;de;; should be stationary for 


v 


given de;; and varying o;; subject to f(o,;) = constant, is 
i 4 i ates 
{a;,;de;;—f dA} = 9, 


where dA is an undetermined multiplier. This is simply equation (6). It 
should be noticed that the requirement of the existence of a true stationary 
value will, for some yield criteria, place restrictions on possible strains. 
In particular it can be shown that the volume change is zero if, and only 
if, the yielding does not depend on hydrostatic pressure, so that f is a 
function only of o;;. Thus the plastic potential is not only consistent with 
the two observed facts of zero plastic volume change, and yielding inde- 
pendent of superposed hydrostatic pressure, but provides a link between 
them. 

[t is interesting to examine the plastic potential idea in relation to the 
behaviour of mild steel. Annealed mild steel yields in accordance with 
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the maximum shear stress criterion of yielding. The yielding is, moreover, 
non-uniform and consists of localized simple shears in the Liiders’ lines 
which traverse the specimen in the yield-point extension range of strain. 
This is exactly the type of deformation corresponding to a potential 
f(o;;) = Tmax; that is, a shear in the direction of the maximum shear stress. 
On the other hand, cold-worked mild steel yields uniformly and obeys the 
Huber-Mises and Lévy-Mises laws which, as we have seen, are also con- 
sistent with the existence of a potential. 


3. The general principle of maximum plastic work 

In the last section a maximum work principle was established for an 
element in free flow obeying the laws (1) and (4). A general statement 
of this principle will now be proved for a finite mass of plastic material, 
no element of which is being unloaded. 

Suppose a plastic mass is in quasi-static equilibrium under a stress 
system o;;. Let u; be the flow velocities on the bounding surface of the 
mass. Velocity in this context is not measured on a time-scale, since 
indefinitely slow flow is being considered, but relatively to the prescribed 
jj Satisfies 
the equilibrium equations, the yield criterion (1), and the free-flow rela- 
tions (4). Then the rate at which the surface forces of this system do 


quasi-static motion of some part of the system. The system o 


work is greater than that for any other system 6,;, for the same surface 
velocities, provided 6,,; satisfies the equilibrium equations and the yield 
criterion. The proof is as follows. 


The actual rate of work of the external forces 
- | o,;u,l,d8 (/; = direction cosines of the outward normal) 


c - 
=: | -(o;;u;)dV (through the volume) 
OX; 
OU; 477 ‘ ‘ ee . ; 
o,;-— a) (in view of the equilibrium equations) 
0x; 
j 


[ Aoi oi; dV 


by using equation (4) expressed in terms of velocities, with A as the new 
proportionality factor. The rate of work for the ¢;; system with the same 
surface velocities is 
_ . Qo. , _. 
[ 4 u,l,dS = { Fy, (is u,)dV = | 6;;— dl 
o ~s 


Ox 3 


(since the system 6;; is assumed to be in equilibrium) 


=f Adiso%, dV. 








Now 


This 
vectc 
the s 
to sy 
Th 
theor 
is eit 
mass 
static 
tion 
veloc 
(9;;5 2 


over 


e 


| ( 


. 


and 


| ( 


. 


Since 


Hen 
and 
dete 
the 

of t 
exte 
pon 








MAXIMUM PLASTIC WORK IN CLASSICAL PLASTICITY 





Now A is essentially positive, and so the principle is proved if 


04,04; > Gyo; Whenever ojo; = Gj 64;- 
| This is equivalent to the statement that the scalar product of the stress 
} vectors to any two points on the yield circle (Fig. 2) is never greater than 
the square of the radius, and is clearly true. The equal sign corresponds 
' to systems differing only by a hydrostatic pressure or tension. 
This variational principle can be adapted to prove certain uniqueness 
theorems. If arbitrary velocities are taken over the surface, then there 
| is either no possible corresponding stress distribution for which the whole 
| mass is plastic, or, if one exists, it is unique (apart from a uniform hydro- 
| static pressure). The first part of this statement follows from a considera- 
tion of special cases; a general criterion for determining possible surface 


} velocity distributions has not been found. To prove the second, suppose 
(c;;,u;,A) and (&;;,i%;,A) could be two possible systems such that u; = a; 
' over the surface. Then 
: | (6,;—9,;)u,1,d8 | Ao;;(G;;—ajj) dV } | N(G;—9;;)(Fj—o;) AV, 
and 
| (o,,—6,,)a,1,dS = | do},(0;;—6;;) dV -4 | AG j;—o%};)(Gj,—o',) AV. 
Since u; = a; on the surface, we obtain, by addition of these equations, 
| (A+A)(6);—9};)(G—o4) dV = 0. 


But (A-+-A) is everywhere positive, and (é;;—0;;)(6;;—oj;) > 0 for all stress 
systems. Hence:a;; = o;;. The theorem can be slightly extended to the 
case when the velocities are prescribed only over a part of the surface; 
over the rest certain external components of force are zero, while those 
velocity components for which the non-zero force components do work 
are also prescribed. 

Similarly, if the stresses are specified over the surface, then either no 
completely plastic state of stress exists, or, if one exists, it is unique. 
For, if o;;, ¢;; could be two possible states, it has been seen that 


( (6,;—o;;)u,l,d8 4 | \(Gj,;—%;) (Gj —94;) AV. 


Hence, if 1;é,; = l,;o;; on the surface, both sides of the equation are zero, 
and the proof goes as before. The velocities are not, of course, uniquely 
determined. A more general statement can be proved for the case when 
the external forces are prescribed only over a part of the surface, some 
of the velocity components being zero over the rest; in this part the 
external force components which do work for the non-zero velocity com- 


ponents are also given. 
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It is interesting to consider the maximum work principle in relation 
to a familiar energy theorem in elasticity, which is invariably stated in 
a minimal form. It will now be shown that the two variational principles 
are special cases of the same, more general, theorem. Suppose two sets 
of quantities o,; and uw, are related by the system of equations 


l/ou, . ou; 
2 y] = , ' 
sla, Ta) = Wy + Bod;;, 


“ Ox; OX; 
where «, 8 may be variable but are positive, and o = 30,;. Again g;; 
Be a ee ; 00%; ; 
satisfies the equilibrium equations —“ = 0. Define 
OX; 
f=? | (xo; 0; +380) dV; J | o,;U,1,d8, 


where / is an integral through a volume, and J is an integral over the 

surface of the volume. Then, if &;; satisfies the equilibrium equations, 
I—J > I—J, 

where J, J denote J, J with é;; substituted for o,; while «, 8 and the u,; 

are unaltered. This statement is proved by observing that 


I—I = | (€;—o;;)(a04;+Bo8,;) dV + 1(6;,;—0;;). 


ie) 
in which /(¢;;—o,;) denotes that ¢,;,—o,; is substituted for o,; in I(o;;). 
Therefore , ; 
I-I= | i dV41(6 
_ G;;—0;;) — aV + (¢;;—9;;) 
J Ox; 


= o— 7 +I (6;; —0;;); 

after using the equilibrium equations and Green’s theorem. But since 
x, B are positive, /(6;;—o;;) is also positive for all o;;, é;;, and the result 
. , ] 1—2v : , 
follows. By setting « 50 p= ~—— and interpreting wu; as the total 
displacement we recover the familiar energy theorem in elasticity. This 
is usually stated in the form: if the boundary displacements are given, then 
of all statically possible stress distributions the actual one makes I—J a 
minimum. In the actual solution J = $J = potential elastic energy. If, 
on the other hand, we set « = A, 8 = 0, and interpret the w; as velocities, 
we recover the maximum work principle in plasticity. For J = J by the 
yield criterion, and so (—J) <(—J) or J >J. There is thus nothing 
peculiar in having a maximal principle in plasticity and a minimal 
principle in elasticity. 

All the above theorems remain true for a material which work-hardens, 
provided elastic strain-increments are neglected. The value of k (or Y) 
is then a function of the previous strain-history, and varies from point 
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to point of the material. In the statement of the variational principle 
the stress system ¢;; must be such that 6j; 6; is not constant but is equal 
to the given value of oj; 0;; at each point at the moment under considera- 
tion. Free flow is not possible for a material which hardens, since the 
stress always changes during the prescribed strain-increment, and the 
elastic strain-increment is not zero. When hardening occurs, therefore, 
the theorems are strictly true only for a fictitious plastic-rigid material. 
4. Application of the variational principle 

The variational principle of maximum plastic work can be used to solve 
certain special problems. However, its range of application is probably 
less wide than that for the analogous theorem of elasticity, since the 
boundary conditions cannot be prescribed arbitrarily if the region under 
consideration is to be entirely plastic. The principle will be applied here 
to determine the stress distribution in a prismatic bar of arbitrary uniform 
cross-section, plastically deformed by certain forces and couples applied 
to the ends. The cylindrical surface of the bar is supposed to be stress- 
free, and elastic strains are neglected. Cartesian axes of reference are 
taken so that the z-axis is parallel to the generators; the axes of x and y 
lie in the plane of a cross-section. For present purposes it is not necessary 
to select any special point as the origin of coordinates. 

[t is first necessary to examine what surface velocities can be prescribed 
if the bar is to be completely plastic. For simplicity consider only condi- 
tions where the stresses and the rates of strain do not depend on z. This 
rules out the possibility of flexural forces, but allows certain combinations 
of bending couple, twisting couple, and longitudinal force. Let wu, v, w 
be the component velocities. It is easy to show that, if the material is 
incompressible and the rates of strain are independent of z, the most 
general expressions for the velocities (apart from a rigid-body motion) are 

“= - :.. 4A (x?+-27)—3C2x+ Dyz; 
oy 
op 


Ca 


4 B(y?+-27)—4Cy— Daz; 


w = &(x, y) + Axz+ Byz+ Cz. 

¢, % are arbitrary functions of x, y; % determines the warping of a cross- 
section. A, B, C, D are dimensional constants related to the rates of 
bending, extension, and torsion. Of these possible velocity distributions 
only a special sub-set represents the deformation of a completely plastic 
bar; this will be examined subsequently. 

These velocities will now be used in the variational principle. In addition 
to the assumption that the stresses do not depend on 2z, it is supposed 
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that o,,, %,, Gz, are identically zero, by analogy with the corresponding 
problem in elasticity. Such assumptions are always necessary to obtain 
a tractable variational problem; their validity must be verified after- 
wards by substitution in the equations (1) and (4). In view of the equi- 
librium equations we can introduce a stress function f(x,y) such that 

22 — eZ, oy; = ra, (8) 


hoy - Ox 


Co 


where f = 0 on the contour of the section since there are no external 
forces on the cylindrical surface. The longitudinal stress o,, is given by 
the yield condition: 

o,, = V3k(1—f2—f?)}, (9) 
where f,, f, are written for éf/éx, éf/éy. The rate of work of the external 
forces is then az 

W | | (o,,u+o,,v+o,,w) dady, (10) 


Syz 

where the integral is taken over the ends of the bar. This expression has 
to be maximized for surface values of uw, v, w (tentatively assumed to be 
given by equations (7)) with respect to stresses satisfying equations (8) 
and (9). It is not necessary to know in advance the actual velocities on 
the cylindrical surface since the external forces there are zero, and con- 
tribute nothing to the work done. 

In equation (10) terms not involving z cancel out at the opposite ends 
of the bar. If L is the length of the bar, the remaining terms lead to 


W att | 
li | | (Av+ By+ C)(1—f2—f})* dady—D ik (af,,t+-yf,,) dady+ 


+4L || (Af,—Bf,) dady. 
On the right-hand side the third term is zero since f = 0 on the contour, 
and the second term can be transformed to give 
W 


LL V3 ia) (Ax+ By+ C)(1—f2—f?)! dady+2D ia) fdady. (11) 


Stationary values of W then correspond to functions f satisfying 
shee see] i‘ ze eae. id 2D _ 9, (12) 
dz /(1—fz—S5) oy (1—f2—f?) NB 
The solution of this equation, subject to the boundary condition f = 0, 
determines the stress function f and thereby the stresses o,., o,,, 9-,- 


YZ 
The equation is elliptic, as can be seen by writing it in the form 


(Ax+ By +O) SD) fn fx Syfayt +(1—f; 2 Soyt + 


. (I 
+(Af,+ Bf,)(1 —fi—fz)+ / —f2—f2)' = Q. (13) 
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The particular integral 1—f2—f? = 0 corresponds to the case of pure 
torsion. Generally it will probably be more convenient to apply Rayleigh- 
Ritz methods in the variational problem of equation (11), rather than 
solve equation (13) directly. 

We have now to verify that a corresponding velocity solution can be 
found such that the Lévy-Mises equations are satisfied. By substituting 
from equations (7), (8), and (9) in equation (4), it is readily proved that 





d(x, y) = 4(Ay®— Ba3)+hay(Axr— By); (14) 
| 3(Av+ By+C)f, 
op Dy —3{42t Byt My, (15) 
Cx V(1—fe—Sy) 
eb _ p, , V3(Aet a. Ir (16) 
Cy \ (1 —fz —fi) 


The condition for the compatibility of equations (15) and (16) for the 
warping function % is just equation (12). Substituting the expression for 
¢ into equation (7), the velocities for a completely plastic state are 


u 1A (y?—x?— 22?)—} Bry—}C2x+ Dyz; 
v } Avy+i B(2?—y?—22?)—4Cy—Daz; (17) 


w = &(x, y)+ Aaz+ Byz+ Cz. 
The corresponding external forces applied at the ends of the bar are 
statically equivalent to a longitudinal force 


Z = 3k {{ (1—f2—f2)' dedy, 
a twisting couple | G. = 25 (  f dady, 
and a bending couple with components 
G, = v3k [ ( y(1—f2—f?)* dady; G, = V3k [fx —f2—f?)* dady. 


It should be noted that the sign of the radical is not necessarily the same 
over the whole cross-section. To given values of the three ratios A/D, 
B/D, C/D, expressing the relative rates of bending, twisting, and extension, 
correspond values of Z, G,, G,, G,. On the other hand, three ratios of 
the latter cannot be arbitrarily prescribed if the whole bar is to be plastic; 
for example, a single bending couple (@,, G,,) does not generally produce 
a completely plastic state. 

Finally let us consider the present variational principle of maximum 
plastic work in relation to the heuristic maximum ‘effort’ principle pro- 
posed by Sadowsky (15) for the ideally plastic body. Sadowsky suggested 
that ‘among all statically determined possible stress dist:;vutions (satis- 
fying all three equations of equilibrium, the condition of plasticity, and 
boundary conditions) the actual stress distribution in plastic flow requires 
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a maximum value of the external effort necessary to maintain the flow’. 
Sadowsky verified this tentative principle in the case of a uniform bar of 
circular section, plastically deformed in combined torsion and tension, by 
maximizing the longitudinal load for a given torque or vice versa. Prager 
(11) showed, more generally, that the principle could be applied to obtain 
the correct differential equation for a uniform prismatic bar of arbitrary 
section under combined torsion and tension. Handelman (3) verified 
Sadowsky’s principle for a bar under combined bending and torsion, the 
section having an axis of symmetry, by finding stationary values of the 
bending moment for a given torque. The reason for the success of 
Sadowsky’s heuristic principle in all these cases (and in certain others) 
is that the actual velocities turn out to be such that the problems of 
finding stationary values of the work done, and of the ‘effort’ (i.e. load, 
bending moment, torque), are identical. This identity has already been 
observed in the special case of combined torsion and tension of a circular 
cylinder (Waters, 20). Generally, of course, Sadowsky’s principle does 
not lead to correct results, even where the applied forces are such that 
an unambiguous meaning can be attached to the term ‘effort’. 
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SUMMARY 

The function ys(k)—log k+1/2k is real for real positive k, and complex for pure 
imaginary k, its imaginary part being then given by a simple formula. Its real 
part, considered as a function of 1/k?, varies smoothly through 1/k? = 0, and is 
a convenient auxiliary function to use for interpolation purposes ; it is also required 
in connexion with the tabulation of the confluent hypergeometric function. 

The real part of %(k)—log k+1/2k has been evaluated, and is tabulated, as a 
function of 1/k? for the range 1/k? 1-00(0-01)+- 1-00, to 8 decimals. 


1. Introduction 

In one group of applications of wave mechanics to the calculation of the 
structure and properties of many-electron atoms, it is required to have 
tables of two independent solutions of the equation 


dy [1 ha |» 


da2° |x ni x 





= 0 (1) 


for integral values of 1, y being a constant which may take both positive 
and negative values, including zero (1, 2). There is a recurrence relation 
between solutions for three successive integral values of | (2), so that it is 
sufficient to consider the cases / = 0 and 1 = 

Solutions of equation (1) can be expressed in terms of the confluent 
hypergeometric function W,.,,(z) of Whittaker (see ref. 3, ch. 16); two 
solutions are 

y Wye144(%/k) and y= W_jers4(—2 k), (2) 

where k? = 1/y. But this form is inconvenient in this context, first, because 
here the variation of the solution with y at constant x (not at constant 
yx) is important, and this is obscured by writing the solution in the form 
(2) (see ref. 1); secondly, the solutions for different values of y are required 
as functions of x, not as functions of y!x, and thirdly, because the form (2) 
is clearly not appropriate for y = 0, and looks likely to be awkward for 
small positive and for negative values of y, all of which are important in 
this application. 

Consider ‘ation of the convenient solutions of (1) to take in this context 
(see ref. 1) leads to the conclusion that in the expansion of one of these 
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solutions near the origin a coefficient has to be used which involves y in 
the combination b(1/yt)+4log y+} yt, (3 
where u(z) is the logarithmic derivative of the ['-function, 
%(z) = I’ (z)/T(z) 
(the notation used here is that of Whittaker and Watson (3), chapter 12); 
for positive values of y, the positive value of y* is to be understood. It is 
convenient to write k? = 1/y (4) 
with the same understanding Then the quantity (3) is, 
wk(k) —log k+1/2k. (5) 
The values of the function (5) are required, as a function of 1/k?, for real 
positive and negative values of this argument. 

For y positive (k real) the function (5) is real; for y negative (k pure 
imaginary), it is complex. From the asymptotic formula for the ['-func- 
tion, which holds provided |jargz| < 47, it follows that 

20 
ys(z)—log z+-1/22 ~ ¥ (—)B,/2rz*” 
r=] 
(the notation for Bernoulli’s numbers here is that of Whittaker and Watson 
(3), $7.2), that is, 
ys(z)—log z+1/22 ~~ — sla mhiae es i x0) (6) 
I2]2* 102% 2126 2028 § 11 210 
so that for small positive values of y, ¥(k)—log k+1/2k behaves like — jy, 
and is zero for y = 0. 

Formally, the asymptotic series (6) does not apply to pure imaginary 
values of z (y negative); but it is found that both the real and imaginary 
parts of %(k)—log k+-1/2k, regarded as functions of y, difference smoothly 
through y = 0. This function is thus a convenient auxiliary function to 
use for interpolating (kh) for real positive values and pure imaginary 
values of k, and for this reason may have applications beyond those for 
which it was originally required; it has therefore been calculated and 
tabulated to a greater accuracy than that required in that context. 

2. y>0 

For y > 0, k = y~* is real, and there is then no difficulty in evaluating 
w(k)—log k+-1/2k, either from Davis’s tables (4) of %(2) (written V(x) in 
these tables), or, for sufficiently small y, from the asymptotic series. To 
obtain, from Davis’s tables, values of ¥(k)—log k+-1/2k for integer values 
of 100y = 100/k?, it is most convenient to evaluate %(a)—log a-+-1/2a for 
a small group of values of x, in the neighbourhood of k = 10/(100y), for 
which (x) can be taken from Davis’s tables without interpolation, and 
which bracket the required value of k, and to interpolate between these 
values, rather than interpolating between the values of x(x) itself. 
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3. y <0 
For y < 0,k yt is pure imaginary, whereas the asymptotic series (6) 
is only valid for |argz| < 47. For pure imaginary values of z, say z = ty, 


[(iy)| = (2/ysinh zy)! 
(Whittaker and Watson, ref. 3, ch. 12, Example 7), and hence the 
imaginary part of (ty) is 
Im %(iy) = $a coth ry+ 1/2y (7) 
(see ref. 5, p. 203), so that for k = i\y- 
Im] (k)—log k+1/2k] = a/(e?**—1) 
m/[exp(27\y-*|)—1]  (y <9). (8) 
This is a function for which all derivatives of finite order are zero at y = 0, 
so that it joins smoothly at y 0 on to 
Im[¢(k)—logk+-1/2k] = 0 (y > 0). 

Thus numerical methods are only necessary to evaluate the real part 
of &(k)—log k+-1/2k for (1/k2) < 0. The method used was to evaluate the 
real part of u(n+i\k|) for n 10, for each negative value of 1/k?, from 
the asymptotic series (6), and then to step down to n = 0 using the 
recurrence relation b(z) = (z+1)—1/z. (9) 
The required tabular value is 

Ref (i \k|)—log ¢\k|+-1/2¢|k|] = Rep(t|k|)—log |k|. (10) 
This was done for each integer value of —50y from 0 to 50. 

If the terms up to that in 1/z!° in the asymptotic formula for %(n-+-2|k]) 
are taken, the error in the result is less than 1.10-" for n = 10, so that 
use of these terms is adequate to give 10-place accuracy. The real parts 


of these terms give 


; n?t+|k\? 1 n 
Re ys(n-+i|k|)—log |k| = }log —L+ + —-= ——__ 
. k|? 2 n?+- |k|? 
lf n?—\|k/? 1 n4—6n?\k\?+ |k\4 
12} (n?+|k|?)? 10 (n?+ |k\?)4 


(the coefficients of the terms in the numerators of the fractions inside the 
square bracket are alternate binomial coefficients). With » = 10 and 


1/k?, this gives 


O}4 
Re w(10+i\k|)—logk = $log(1+100}y =. =. —_ 
' 2 1+-100\y 
l 100\y|—1 1 ,(100}y|)?—6(100\y|)+-1 
~ yo 4 (100|y|+-1)? cs i ~ (100/y|-+1)4 ond 
l , (100!yv})8?—15(100 21. 15(100\y|)—1 
2) ee ee. (11) 
21 (100|y|+1)8 


and the real part of (9) gives 
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Re| uk(k)—log k+-1 2k| 








R ] 2 3 9 
Xe Y(10-+-7|k|)—log |k -lisy ie ita r atari Ti 2 | 


; ] 2 9 
= Re %(10+7 k )—log k aaah Fes ate <5 Be ii (12) 
— | e 17 T 


Since the calculations were done for exact values of y, the forms (11) and 
(12) were very convenient for computation. 


4. Results 
Values of u&(k)—logk+1/2k were evaluated to 10 decimals for 
y = 1/k? = —1-00(0-02)+1-00, and differenced. The differences indi- 


cated that random errors due to rounding-off were not more than 2 in 
the tenth decimal, although, for negative values of y, each value was 
calculated as the sum of about 12 terms. The values were subtabulated 
to 0-01 intervals in y and rounded-off to 8 decimals; the results were 
differenced on a ‘National’ machine and are given in the table. The 
maximum error in a tabulated value should not be more than 0-6 in the 
eighth decimal. 

The first and second differences of the function values are given in the 
table; it is never necessary to take fourth differences into account, and 
linear interpolation is adequate to give six-decimal accuracy. 
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FINITE DIFFERENCE FORMULAE FOR THE 
SQUARE LATTICE 


By W. G. BICKLEY (Imperial College, London) 
[Received 20 October 1947] 


SUMMARY 

The paper gives approximate formulae for derivatives (including combinations 
like V? and V*), and integrals, of a function of two independent variables, in terms 
of its values at nodes of a square lattice, primarily for use in the numerical solution 
of partial differential equations. Consideration is given to the form, as well as to 
the magnitude, of the leading terms in the error, and what is believed to be for 
most purposes optimum combinations are thus selected for the simpler compact 
sets of nodes. 


1. Introduction 


THE use of numerical methods in solving problems in applied mathematics 
is becoming increasingly common for partial, as well as for ordinary, 
differential equations. The first step in such a process is usually the 
replacement of differential equations by their finite difference approxima- 
tions. For ordinary differential equations, approximations of various 
orders of accuracy to derivatives and integrals are well known and listed,+ 
but for partial differential equations it is as yet rare to find any but the 
crudest approximations employed. {§ 

For many purposes formulae of the Lagrangian type are convenient. 
In dealing with partial differential equations with two independent 
variables, equidistant, and equal, intervals in the independent variables 
will usually be employed, so that the required function is computed at the 
nodes of a square lattice (net). The object of this paper is to explore the 
possibilities of refining approximations to the quantities—derivatives and 
integrals—which occur in the formulation of problems in applied mathe- 
matics. 

Considerable use is made of symbolic operators, which are powerful in 
developing such formulae; a posteriori proofs of the formulae are readily 
constructed—if they are deemed necessary. 


+ W. G. Bickley, ‘Formulae for numerical integration’, Math. Gaz. 23 (1939), 352-9; 
l., ‘Formulae for numerical differentiation’, ibid. 25 (1941), 19-26; L. J. Comrie, Inter- 
polation and Allied Tables, H.M. Stationery Office (reprinted 1942). 

t Some attempt is made by A. Thom, ‘Arithmetical solution of equations of the type 
V‘t = const.’, A.R.C., R. and M. (1933), 1604. 

§ Formulae of various types are given by L. Collatz, Ligenwertproblem und ihre numerische 
Behandlung, Akad. Verlag. (Leipzig, 1945). 
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We refrain from giving examples of the use of the formulae since this 
is quite straightforward. Nor do we advocate that the refinements should 
be used in the early stages, or on the coarser nets. Their proper place is 
in the final stages, where they may add considerably to the accuracy of 
a solution with but little extra labour and without further reducing the 
mesh length. 


2. Symmetric sums 
We cover the (x, y)-plane with a square lattice of mesh length a, and 
label a typical point 0 (see Fig. 1). Neighbouring points are numbered 
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i 12, as in the figure, and these numbers will be used as subscripts 


to indicate that the value of any quantity (usually f(a, y)) is to be taken 
at the point in question. 
We commence by recalling the symbolic form of Taylor’s series in one 
variable, 
f(ath) = e*P f(x), (1) 
where D = d/dx. 


To deal with two independent variables (x,y) with a mesh length a it 
is convenient to introduce the symbolic operators 
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In terms of these we may write 
hi ef, te = enfy ) 
h=ef. f=ee | (3) 
Is es "fo; 


and so on. 


We note that €2+ 7? = a®V? (4) 
and find it convenient to write 
; , & adi ; 
o7 = oe — = @*9?*, (5) 
OXCY 


Now it will usually happen that the quantities occurring in the mathe- 
matical equations will be independent of the choice of axes—invariant 
for rotation of these axes—and clearly values of such quantities at 0 can 
be represented or approximated only by sums of values at points sym- 
metrically disposed about 0. Of such ‘symmetric sums’ three only (apart 
from f, itself) are simple (in the sense that they include only four terms), 
and involve points not too remote from 0, namely, 


S,=fitletSsths (ef+e7+e-$+¢e-") ie= = 2(cosh €+ cosh ») fy 
4fot (€+ y*) fo tF “(+a fero (€6+- 9°) fot+ 
4f,+a°* "V2 ‘fot of seesalh uitihtialitall, . 


(6) 
So iP Lf T Ja Is 
(ef +n | e—§+nt e-§-n 1 1) fo 
4 cosh € cosh 7 fy 
; a a ee 
4fo+ —— n”)to+ .(§* + 6F"n* + *) )fo+ 
Z 4! 
| 4 (€8 4 r£2.,4_! 6 ! 
+ (69+ 15649? + 1569+ 9°) fot 
> ahh + ta4*(V4+-4D4) fot qiqa®(VS+ 12V?D) fot... - (7) 
Ss = fothiotSuthe (e% + e7+-¢ 26 +-e-2) fy 
sf, i 10°V2f, T sa4(V* 2 *) fo T is(V8—3V?2D4) fy+.... (8) 


The next symmetrical sum involves eight terms: 
f(+2a,+a) and f(a, +2a). 
3. Approximations to V2 
By equation (6) we have immediately the well-known and most fre- 
quently applied approximation to V?f, namely, 
of Lig? 2 
V*fo (S,—4fo}/a?+O(a?), (9) 


usually attributed to Liebmann. 
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Equations (7) and (8) lead to what is essentially the same result for 
mesh lengths av2 and 2a respectively. 

The problem of improving the accuracy conveniently is not simple. 
Clearly the use of S, and S, cannot eliminate both V4 and @*, and although 
we can eliminate both by using S, and S,, we use points distant 2a from 0, 
which complicates procedure at a boundary, and also we cannot then 
obtain any help from §,. 

The use we can best make of S, and S, is to eliminate the term in 9, 
since this is not invariant for rotation of axes. The resulting ‘error’ term 
in a* then involves only V‘f,, and is thus (a) invariant, (b) calculable (at 
least approximately) from the values of V?f, and (c) small (theoretically 
zero) for Laplace’s equation. 

The resulting formula is 

48, +S, = 20fo+ 6a°V fot daV “fot goU(VI+ 2D "fot... (10) 
and from this} 
V2fy = {48,+S,—20f,} /6a2—3,0°V4f,+O(a*). (11) 

One may use (9) to compute V‘f, from the (approximate) values of 
V?f, with an error of order a?, so that if the term —;,a°V‘f, is allowed for 
in this way, the error in (11) is of order a*. This procedure is equivalent 
to the use of a fourth difference correction (as advocated by Fox),{ and 
has, indeed, the advantage that it is a two-dimensional correction, whereas 
all Fox’s fourth (and higher) differences are differences of one-dimensional 
sequences. 

We can, alternatively, eliminate completely the a* terms between (6) 
and (8), obtaining 

16S,—S, = 60f)+ 12a°V2f,—,a9(V4*—3F4*)V2fy... (12) 
or V2fo = {16S,—S,— 60f,}/12a?+O(a!). (13) 

If f, is an approximate solution of Laplace’s equation, then the last 
term written in (12) shows that the ‘error’ of (13) is of order higher than 
a*, But, as already indicated, the superior accuracy of (13) is purchased 
by wider spread of the points used. 

Using S,, S,, and S, we may eliminate V*/, and Z*f,, and obtain the 
well-known approximation, 

V*fo = {20fo—8S,+ 28,+83} a*+-O(a?), (14) 
and clearly nothing better is available unless we include additional points. 

+ This formula is given, effectively, by P. M. and A. M. Woodward, ‘Four-figure tables 
of the Airy function in the complex plane’, Phil. Mag. (7) 37 (1946), 259. 

t L. Fox, in his paper ‘Some improvements in the use of relaxation methods for the 


solution of ordinary and partial differential equations’, Proc. Roy. Soc. A, 190 (1947), 31-59, 
corrects the crude finite difference approximation by the use of higher differences. 
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4, Other derivatives 


First-order derivatives occur infrequently in the governing differential 
equations, but are common in boundary conditions. 

The crudest approximation to (@f/éx), is Euler’s, namely, (f,;—fo)/a, but 
this is an approximation ‘centred’ at the mid-point of the mesh 01. For 


the value centred at 0 we use 


fits 2 sinh € fo 


el Ri cies ” 
2é+ oe tee t So (15) 
ov: 
so that (of /Cx), (f:—fs)/2a+-O(a?). (16) 


To improve this we must clearly use approximations to (éf/éx), and 
] : 2 
(ef/éx),, in equal proportions. Now 
fs—fetfa—hr = (e&+7—e-$+0+ e& 1-4-1) fy 
4sinh € cosh nf, 


IF Tas tT - 
Hl Fetoitt-|efe (17) 


It is clear that potentially the most useful combination will be one in 
which.(since they cannot be eliminated) the third-order terms constitute 
a multiple of V2(éf/éx). We readily find the combination 

4(f; ts) 1 (fs Ts) I (fs J7) p12 T 2(€*+- y?)+..j}Efo 
1 2a(ef/éx) 9+ 2a°V?(ef/ex)y+O(a*).... (18) 

By using the approximation similar to (9) in the second term on the 
right-hand side of (18), the relative error in (éf/@x), may be made of 
order a‘. 

When the point 0 is a boundary point, the use of the above formulae 
necessitates values at ‘fictitious’ points, outside the boundary, but the 
use of such points is a commonplace method of the ‘relaxation’ and other 
techniques. 

Soundary derivatives at points other than nodes, and normal or tan- 
gential derivatives for curved boundaries, do not lend themselves to any 
simple treatment if any improvement upon the crudest approximation is 
desired. 

The formula for (éf/éy)), corresponding to (18), may immediately be 


written down. 
5. Surface integrals 


In many problems the integral of a function taken over some region, 
determined at the nodes of a lattice, is required. 
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The simplest approximation (corresponding to the trapezoidal rule in 
one dimension) associates the value of the integrand at the node with the 
square of side a centred at the node (Fig. 2), and applies a factor } to 
nodes on a boundary and } to nodes at corners. 
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An obvious improvement is to use Simpson’s rule in both directions, 


According to this 


a a 


, o£, = — 
| | fdady = 5 (16fo+ 48, +53), (19) 
ani me ; 
but this is not the best that can be done with the nine nodes 0, 1.,..., 8. 
- We have 
s 3 , 3 
| | f(x,y)dady = a* | | er$+snf drds 
a “a = . 
»sinh é sinh 7 , 
= 40° g 7 


En J0 
9 ' 2 2\ 1 l 2,2 4 
= 4a? 1 Stn) (E+ Pent nt) tl 


= 4a*{1+ fa°V2+ 54,04(V4+$D4)-+...} fo. (20) 

By the use of f, S,;, and S, we can obtain the correct amounts of any 

three of fy, V2fy, V4fp, and Z*f,—but not of all four. If one must remain, 
clearly it should best be V‘4/,. We find, with no difficulty, 


a a 


[ [ f(e,y)dedy 75 (88fo+ 165, + 752) —a°V4fy+ O(a’). (21) 
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The numerical coefficients have become larger—beyond the scope of 
mental arithmetic—but the formula has the merit that the value of the 
leading term of the error is readily calculable, and is small for solutions 
of the Laplace equation. 

We may also compare (19) with (20), and the result is 

a . . 
[ f f(x,y) dandy = “ (16f,+48,+S,) (V4 294) f,+-O(a*), (22) 
eS 9 45 

in which the term in a® is not readily computed. 

Two cruder formulae may be mentioned, giving approximations to the 
same double integral, namely 


3a>(2f,+S,)+O0(a'), (23) 
4a*(8f,+S,)+O0(a*). (24) 


Finally, it may be worth recording the approximation to the integral 
over the square of side a centred at the point 0 (Fig. 2), in terms of the 


values at the points 0, 1.,..., 8. 
ha 7 
a | , sin 1d €sinh $7 
| J (x, y) dady 42 ; If, 
ye $7) 


a*\1+3,a* my" + - 70506 u4(Vi4 $F") 9 So» (25) 


and comparison with (6) and (7) leads to the result 


( [ fw. y) dady 1244f,+ 38S, + 118,)—3i,a°V4f,+O(a8 


{ 
| 1440 
la 4a 1 (26) 


It is at first sight attractive to apply this successively to all nodes of 
a region. The result is that the factor for all internal nodes is unity— 
which clearly cannot be expected to give the best results. The explana- 
tion lies in the fact that the regions to be associated with nodes on the 
boundary are not complete squares but half-squares or (at corners) quarter- 
squares. The contributions from these are not given by symmetric sums, 
so that the corrections must be applied to a border of nodes on and within 
the boundary. The situation is parallel to the Gregory integration formula 
in one dimension. Corresponding correction terms could be worked out, 
but their use would be a complication which (21) or (22) completely 


avoids. 
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6. Computational molecules 

It seems serviceable to represent pictorially the formula given above, 
in the form of diagrams indicating the factor to be applied to the function 
value at the corresponding node. In each case the result of adding the 
given multiples of the function at the corresponding points is given, but 


for the ‘error’ terms reference must be made to the appropriate equation. 
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ON LAMINAR BOUNDARY-LAYER FLOW NEAR 
A POSITION OF SEPARATION 


By 8. GOLDSTEIN 
(Department of Mathematics, The University, Manchester) 


[Received 27 October 1947] 


SUMMARY 

Singularities are considered in the solution of the laminar boundary-layer equa- 
tions at a position of separation. A singularity of the type here considered occurred 
in a careful numerical computation by Hartree for a linearly decreasing velocity 
distribution outside the boundary layer; it may occur generally. Whenever it does 
occur, the boundary-layer equations cease to be valid at and near separation on the 
upstream side, and also downstream of separation. The work suggests that singu- 
larities may arise in the solution of non-linear parabolic equations due to their 
non-linearity. The formulae found may help computers of laminar boundary layers, 
who desire more than a rough solution, to have an end-point at which to aim. 


1. Introduction and summary 

For a flow at a large Reynolds number along an immersed solid surface 
a boundary layer is formed through which the velocity rises rapidly from 
zero at the surface to its value in the main stream. The approximate 
equations for the two-dimensional flow of a fluid of constant density p 
and kinematic viscosity v in a boundary layer are 


ou ou 1 op | onu 
u -v —— +9 — 
Ox oy p ox cy (1) 
Cus Gus 
u P v — 
oy Ox 


where x is distance measured along the solid boundary in the plane of 
the flow, y is distance normal to the surface, wu, v are the velocity com- 
ponents in the directions of x and y increasing, p the pressure, and & the 
stream function. According to the approximations of boundary-layer 
theory, p and @p/éx may be taken independent of y, and if U is the 
velocity just outside the boundary layer in the main stream, 


l Op pau 


_— 
bo 
— 


p ox dx ; 

Moreover, according to these approximations v/u and (@u/éx)/(@u/éy) are 
small 

As boundary conditions we have that u = 0 and v = 0 (or & = 0) at 

0, u is given as a function of y for some initial value of 2, and the 
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velocity passes over smoothly into the velocity of the main stream, i.e, 
u—> U, eu/ey > 0, eu/dy? > 0, etc., as y >. 

If dU/dx <0 and ép/éx > 0, then du/éy at the wall y = 0 decreases 
as x increases until it vanishes; beyond the section at which it is zero, a 
slow back flow sets in along the wall, and the boundary layer separates 
from the surface. 

Let.d be any representative length of the system, U, a representative 
velocity, such as the undisturbed stream velocity, and R the Reynolds 
number U,d/v. The equations may be made non-dimensional by writing 


a’ =2/d, y' = Ryjd, uw’ =u/U,, v' = R/U, py’ = p/pU3. (3) 
In the non-dimensional form, let x’ = 0 be the initial section, at which ’ 
is given as or approximated by a polynomial or power-series 
, , 19 
wu’ = a,y'+a,y'?+..., (4) 


Op’ , “ ; 
and let — 2 = Pot p, x’ +pox'?+.... (5) 


Ox 
Then it is known that if there is to be a solution without singularities, 
certain equations must be satisfied: . 
é — - mR! L« — 29 . a > 2 
2a.+pPp = 9, Ag=0, 5!a;+2a,p,=—90, 6lag—2pop, = 0, etc. (6) 
Only a@,, 44, @;,... are at our disposal. When the conditions are broken, 
the solution has an algebraic singularity at x = 0 (1).T 
At the position of separation éu/éy = 0 at y= 0, ie. a, = 0. The 
conditions for the absence of singularities when a, = 0 are considerably 
more complicated than those above.{ If we suppose w expanded in a 
power series in x (we drop the dashes for the present) 
U = Up tu, r+U,2?+..., (7) 
where Up, U,, U, are functions of y, expressible as power series, 
Uy = Ag y*+a, y+... 
u, = b,y+bey?+... }, (8) 
Us = Cy +Cgy?+... 
the conditions are 
») e _ _— sant = ! — a 
2ao.+pp = 90, ag=0, ag=90, as=—0, Glag = 2p, | (9) 
a,= 0, ad, = I, etc. 


Only dg, G9, Ag, Ag9,--- are at our disposal. In addition 6,, c,, d,,... are 
determined, not from the equations for u,, uw, Us,... respectively, but from 


+ In the last term in equation (5) on p. 4 of ref. (1), the denominator should be 8a}, not 
4a}; I am indebted to Prof. Hartree for this correction. See also ref. (2) at end. 

t Goldstein, loc. cit., pp. 17, 18. [In the last line of p. 17 of ref. (1), in the equations 
of the footnote, for 6! ag = 2p» pg, read 6! ag = 2p 9 p,.] 
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the conditions for the absence of singularities in w,, us, U,4,... respectively. 
There is also an ambiguity of sign, which can be determined only from 
physical considerations. If the conditions are broken there is a formal 
solution for the flow downstream of the form 


b= Pl foln)+éhila)+--] | 
w= LL fon) +Ehiln)+-] 
where £ == #, n = y/4at. (11) 


(10) 


This formal solution fails, however, in certain circumstances, one of which 


is that the condition yp 0 (12) 
ee Ge! ie ™ 


is satisfied, while the other conditions are not satisfied. 

No other work has been reported on possible singularities at separation. 

No analytical solution is known for a boundary-layer flow involving 
separation, and the methods used are approximate and numerical. The 
published methods of computation are rather rough, but recently more 
exact methods have been suggested and tried. The work described here 
arose out of an unpublished communication from Professor Hartree, in 
which he repeated Dr. Howarth’s computation (3) for a linearly decreasing 
velocity distribution, U = £,—f,x, with u= U at x=0. Professor 
Hartree replaces the partial derivatives with respect to 2 by finite 
differences, and retains the y-derivatives, so the partial differential equa- 
tion is replaced approximately by a sequence of ordinary differential 
equations, each of which relates the velocity distribution through the 
boundary layer at one section to that at another section a short distance 
upstream, where it is known. The ordinary differential equations were 
solved laboriously on hand calculating machines rather than on the 
Differential Analyser in order that more significant figures might be 
retained. 

Now all computations in which any attempt was made to obtain real 
accuracy at and near separation seem to have met with considerable 
difficulty. As a result of his computations, Professor Hartree was con- 
vinced that there was a singularity in the solution at the position of 
separation, and I undertook to try to find some formulae that would 
hold near this singularity and would help in finishing the computation. 

To study the singularity near separation, the equations are put into 
non-dimensional form in a special way. Let z,, U,, U, be the values of 
x, U, dU /dx at separation, so that U, > 0, U", < 0. We are not interested 
in any other properties or dimensions of the system, so as representative 
length 7 and Reynolds number R we take 


ln—UsR, B= (13) 
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We are also concerned with the flow upstream of separation, so for our 
non-dimensional distances we write 


4, = (x,—2z)/l, y, = Ry/l. (14) 
Also put 


U, = U/U,, 4 = R/U, U,= UU, py = p/pU2, = R*p/1U,. 


(15) 
The equations become 
c ‘ c a9 
du eu Op, , mu 
—uy— +9, = 4 
6x, dy, Ox, oy? 
0 Ous , 
« hoon ey < Tea aa . (16) 
oy, 6x, | 
op, Uy. dU, 
8h | 
6x, dx, 


It is easy to see that U, = 1 and dU,/dx, = 1 at x, = 0, so we may 
write 


“Pi — (14 P.a,+P,22+...), (17) 
OX, 
i.€. Py (in our previous notation) = —1. For the linear velocity distribu- 
tion U = B,—f,2, 
(te), R=1, R-R=-=0 (9 
Ox, 


Formulae for the P,, P,,... are easily found in the general case, and if U is 
taken as a given function of 2, it is easily found that the P’s are inde- 
pendent of the position of separation if 
U = (constant)e-®*, or U = (By—B,2)" or (By+P,2)-™, (19) 
the constants B, By, B,, m being positive. (However, since 1 = —U,/U;, 
the scale varies as x, varies in the last two cases.) For other values 
of U, the P’s depend on z,, the position of separation. 
The boundary conditions are , = 0, u, = Oat y, = 0, and u, > U,, etc., 
as y,>00. Since x, = 0 is a position of separation, (8u,/Ay,),,-9 = 0 at 


x, = 0, so 2 | 3 
1 U, = 2,yi+ayi+... at x, = 0. (20) 


Singularities in the solution for the corresponding system of equations for 
the motion downstream have been considered (see equations (10) and 
(11)); near 2, = 0, y, = 0, ¥, is a function of 2} and y,/x}. The skin- 
friction is (éu/@y),-) and the determining quantity is (Ou, /Ey,),,-9 Which 
is an ascending series of powers of 2}, beginning with a multiple of 2}. 
If, however, 2a,+-p) = 0, which corresponds to a, = } in the new nota- 
tion, there are special features in the solution; in particular the series for 
(0u,/Cy,),,-9 now begins with a term in 2}. Now Professor Hartree was 
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quite certain that this particular feature was present in his computed 
solution. The calculations reported here therefore rest on three assump- 
tions, all of which were satisfied in Hartree’s numerical solution: (i) there 
is a singularity at separation; (ii) there is a finite value of u, at separation 
for y, ~ 0; (ili) ag = 4. Related to (iii) Professor Hartree found (em- 
pirically) that in his solution (0u,/é@y,),,-9 behaved near x, = 0 like a 
multiple of a}, where r is certainly less than 1 and greater than }. Thus, 
we must take 


u, = 4y?-+ asy>+... at x, = 0, (21) 
and as a result we find that 
‘ 22 at 3 | ” . oC< 
(OU; /CYy)y,=0 = 2°(% A Xo 1} + ag X, +a, 24+...), (22) 


where the a’s are constants. (The factor 2! is inserted to conform to the 
notation in § 2.) 

The first purpose of the calculations was to find the connexions between 
lz, Ay, As, Ag,... ANA a,, Xy, Xg,... and other formulae for u at and near x, = 0 
to see if the results fitted the numerical values for the special solution. 
There is no mathematical proof that a solution exists with singularities 
of the type considered near separation, but with the above assumptions 
it is difficult to see how the solution could be of a different type. 

We may remark that in assuming that a, = 3, we are in effect assuming 
that (€?w,/@y7),,-9 is continuous at 2, = 0, and then (0”"u,/éy7),,-9 is found 
to be continuous at x, = 0 for n = 1, 2, 3, 4 and discontinuous for n = 5 
and 6 and probably for all n > 5, though 6”, /éy? is continuous for y, ~ 0. 
More important, it is found that at separation v, and @u,/éx, become 
infinite in such a way that a}v, and x} éu,/éx, have finite non-zero limits 
as x, > 0 for all non-zero y,. The basic assumptions of boundary-layer 
theory therefore do not hold at and near separation. Nevertheless, large 
cross-velocities are to be expected at separation, otherwise the assump- 
tions of boundary-layer theory would not break down. 

The formal solution for the motion upstream is found by writing 


ees x4, 7 = y,/2*2, (23) 
by = 286 fo(n) +Efi(n)+€fo(n)+-.-], (24) 
Uy = 28 fo(n) + Efi (n)+€7fe(n)+-.-.] (25) 


in (16), and equating powers of €. Since %, = 0 and u, = 0 at y, = 0, 
f,(0) = f;(0) = 9, and from the value (21) of u, at x, = 0 we find the 
condition ft 

lim -—. = 2¥a,,. (r= 0,1, 2....). (26) 


saa Ml 


The solution for f, must have a double zero at the origin, and must not 
involve exponentially large terms as 7 00. 
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The condition u, > U, as y, > © is satisfied for x, > 0 if it is satisfied 
at x, = 0, ie. if u,; > 1 as y,; > 0 ata, = 0. 

With a, = 3, the solution for f, is found to be f, = 73/6. The solution 
for f, is f, = a, 7”, and 


>. ‘ 
A, = - lim i 0. (27) 
3 9 3 
VG /] 
T ae a 98 
Then Seo = % TR” ; (28) 
vo 
and a, = —}o?. (29) 


|The a’s and a’s are as in (21) and (22).| The equation for f, becomes 
fs —30°fs +in°fs—8nfs = 5fife—Thfet4hfe- (30) 
The equations for all succeeding f’s are non-homogeneous linear equa- 
tions, with the right-hand side rapidly becoming more and more compli- 
cated; thus for /, it involves /,, f,, fs, P,, for f; the first four f’s, and so on. 
The complementary functions involve integrals of confluent hypergeo- 
metric functions; the particular integrals are very involved. The condition 
for the absence of exponentially large terms in f, ist 





Qtort 
Me am 5413 a, (31) 
5(4! 
1 then f at 
and then from (26 a, = — —_—_qi. 32 
( ) 5 40(4!)2 1 ( ) 
The condition for the absence of exponentially large terms in f, is 
3 
7 
Ne = ———- (3 — 8/2) a3, 33 
. 40091) Jos 9) 
, l 77? P 
and from (26) a, = |-—- oo vf — —2 , (34) 
9 600(4!) 360 


The condition for the absence of exponentially large terms in f,; gives a, as 
a multiple of af, though the constant must be found numerically, and then 
a, is found; but the condition for the absence of exponentially large terms 
in f, does not give a,; it requires that 


fo 


10 


ad 6 9 4 
| Hx? —% +7)exr(—] dn = 0, (35) 
0 


where H, is a complicated function of 7, involving f;. 

Again, a, is determined from f, in terms of a,, a;, and P, a, from f, 
in terms of a,, a;, and P,, and so on until we come to f,,. It is possible, 
though it has not been proved, that a; is determined from fo, «, from fi, 


t+ x! is written for [(2+1). 
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and so on. If so, then only a, remains to be determined, and a, (and 
therefore a,) is probably determined by the condition u,—> 1 as y, +00 
at z, = 0. If so the whole solution is determined at separation. In fact, 
if it is true that all the other constants are determinate in terms of a, and 
the P’s, there is a solution only if it is possible to choose a, so that the 
condition wu, > 1 as y, > at x, = 0 is satisfied. Unless this condition 
is satisfied for every value of a,, it will presumably fix a, in terms of the 
P’s. If a, is so fixed, the non-dimensional velocity distribution at separa- 
tion, x, = 0, and just upstream of separation, for small positive values 
of z,, is fixed in terms of the P’s. Suppose now we have a problem in 
which U is a given decreasing function of x, and wu a given function ,of y 
for some 2 < 2,. There are some U’s for which the P’s are independent 
of x,; otherwise they vary with x,. When separation takes place, the non- 
dimensional velocity distribution at and near separation is independent 
of the initial distribution of uw for the former values of U, and for the 
others it is the same for all initial distributions of uw for which separation 
takes place at the same value of x. This suggests that what has been 
found is an asymptotic solution at and near separation, and that the full 
non-dimensional solutions in the above cases all behave asymptotically 
in the same way near separation. 

It appears that the singularity at separation is due to the non-linear 
character of the equations. It is possible to simulate the phenomenon of 
separation by a linear system of equations, and there is then no singularity 
at separation. For example, the solution of 


Cu O7u — 
— mF Y- =, U latt 0, u=Oaty=0, u finite asy>o, 
( oy- 
(36) 
. ll t y om 
8 u Ly? (1—t—1y?)erf-~-—+2 exp{ —‘— (37 
ad aie 7) Mann pas aT) i 
vhere erf x - ee” dw, (38) 
V7T 
0 
80 (Cu/CYy),,—9 0 at ¢=}. (39) 


We should remark also that in equations that are linear but are other- 
vise similar to the equations here considered (the equation for the 
temperature in the theory of the conduction of heat,+ for example), if 
ve attempt to work backwards (i.e. to solve for negative time) from a 


The relationship of the boundary-layer equations to the equation of heat conduction 
s been stressed by Prandtl (2) (loc. cit., pp. 79, 80) in connexion with difficulties to be 


xpected when u < 0 


5092.1 E 
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singularity we encounter exponentially large terms. With given initial 
and boundary conditions, however, the solution for such a linear equation 
is free from singularities for positive non-zero time, whereas the basis of 
the present discussion is the assumption that singularities may occur at 
separation in the solution of the non-linear equations considered. 

The special case considered by Professor Hartree is one in which by 
a correct choice of scale the P’s may be made independent of x,. There 
is always the possibility, therefore, that the occurrence of a singularity 
at separation is restricted to such cases. Another possibility is that a 
singularity will always occur except for certain special pressure variations 
in the neighbourhood of separation, and that, experimentally, whatever 
we may do, the pressure variations near separation will always be such 
that no singularity will occur. 

It is a necessary consequence of the discussion of the motion upstream 
of separation that a, is negative or zero. Professor Hartree finds a negative 
a, from his special numerical solution. When we consider the motion 
downstream of separation in a similar way, we find that when a, is negative 
the solution downstream is not real. When there is a singularity at 
separation there is no real solution at all farther downstream. When 
a, = 0 there is a solution downstream, but then we have a case in which 
the whole solution is free from singularities. These cases include that in 
which (éu/é@y),,9 = 0 for all x. There must, of course, be restrictions on 
the pressure distributions in order that this should happen, and these 
conditions arise from the condition that wu, > 1 as y, > 0 at x, = 0. This 
our method does not permit us to discuss, but one solution (due to Falkner 
and Skan, 4+) is known in which U = cxv™, m = —0-0904 approximately. 

As far as numerical values are concerned, and comparison with the 
computed values, Professor Hartree fitted his solution to the formulae 
here obtained, and considered that he had obtained a reasonable fit. The 
matter has been recently reconsidered by Mr. C. W. Jones, who has 
tabulated f;, f,, f; and has found that within the accuracy of his com- 
putation, the integral condition (35) for the absence of exponentially large 
terms in f, is satisfied. 

Mr. Jones has compared the skin-friction, the velocity distributions at 
separation not far from the wall, the transition to main-stream conditions 
at separation, and the velocity distribution just downstream of separation 
(at (88,/B,)a = 0-956, where U = B,—f, 2 and (88,/B,)a, = 0-959). A satis- 
factory fit is obtained with a, about 0-47 or 0-48. A satisfactory transition 
to main-stream conditions seems to be obtained, but it is not sensitive to 


changes in a. 


+ See also (5). 
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If we assume that Mr. Jones’s numerical work is sufficient to answer 
certain questions, and to make it plausible that our formulae fit the 
solution in the case considered, we still do not know for certain that a;, 
4g... are determined from the equations for fy, f,3,..., and, if they are, that 
a, (and therefore «,) can be determined from the condition u,—> 1 as 
jp > at ay, 0. It is clear that an adequate discussion is not possible 
by the method used here. Three more important questions also remain: 

(1) Is it correct that the formulae represent an asymptotic solution at 

and near separation ? 


What are the most general restrictions on the pressure distribution 


in order that solutions should exist for which (@u/éy),-9 = 0 for 
all x? 
(3) A singularity when U = 8,—f,x being assumed, is the occurrence 


of a singularity restricted to cases in which the P’s are independent 
of x,? Or does a singularity always occur unless (éu/éy),9 = 0 
for all x? Or does a singularity always occur except for certain 
special pressure distributions near separation, and are experimental 
pressure distributions always of the special type? 

It should be remarked that although there is a certain physical plausi- 
bility in the notion that large cross-velocities should occur at separation, 
the existing experimental information is insufficient to settle the question. 

The work described may be summed up by saying that it throws doubt 
n the validity of the boundary-layer equations at and near separation 
on the upstream side, and also downstream of separation; inferences from 
these equations in these regions, which are fairly common in the literature, 
are therefore also in doubt; mathematically the work suggests that singu- 
larities may arise in the solution of non-linear parabolic equations, due to 
their non-linearity; and formulae have been found which may help com- 
puters of laminar boundary layers, who desire more than a rough solution, 
to have an end-point at which to aim. 


2. The solution upstream 
Substitute (23), (24), and (17) into (16), and equate coefficients of powers 
f€. The equation for fj, obtained from the coefficient of £°, is 
So —3fofot+2fo” = 1. (40) 
Since ys, 0 and u, 0 at y, 0, 
f{0) = f,(0) = 0 (¢ = @,1,4....). (41) 
When €+> 0, » +00 if y, 4 0, and since 2!£n = y,, limu, is given by 
21) if f iad 
lim —* PG. = O43. (42) 


7 » r — 
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The condition that the velocity should pass over smoothly into the velocity 
of the main stream will be considered in §3, when the solution for large | and 
values of y,/2} is considered. 








Since a, = 4, the solution for f, is 
fo = 7°/6. (43) 
The equation for f, is then found to be | The : 
” sod l/(y_] 2’ »| § ee: . 
Sf, — 37 fr + 3(r-+4)y ‘fp —(r+3)nf, seis G,, (44) when 
_ As 
c = Y ; oe ” Qf'2 a Fa ts 
where G, = 0, G, = 4f,f/1—3f1?, (45) 
. equat 
and for r > 2, expai 
— , , P 
a * Sy. on Lf a a) ee P | . ; 
G, 2 Le ae 3) fstr 3 (r—2s-++ 2)fsfr-sIt+ Pu, (46) j 
8 
= , ; — - - ff respe 
P.,, being put equal to zero except when }r is integral. The solution for WI 
jf, with a double zero at the origin is 
J, = um 7°, (47) i, 
where a, is a constant; hence from (42), 
Ee ; 
a, = — lim hi = 0. (48) 
Va no 7)" so 
The solution for f, with a double zero at the origin is now found to be Ri 
2 “1 
So = %4?—-— 9, (49) 
0 
and 
where a, is a constant; hence from (42) Ri 
tas 
oe . 
a, = —Fay. (50) 
In order to write down the general solution for f, with a double zero at Henc 


the origin, and to consider its behaviour as 7-00, some discussion is | 
necessary of the complementary functions, and it is advisable to break | 


off and discuss generally the complementary functions of the equation 
a ) and 
for f,. 


Three independent complementary functions are 7*, g,, and h,, where, 


if the function ,F,(a,b,x) is defined by oo 

s aa 19 

F(a,b,2) = 14 @ a a(a+1) silts a(a+1)(a+2) w+... (51) | 
1.b 2!b(b+-1) 3!b(b+1)(6+2) But! 
“ 3__ 1p)! 1! »4m+1 ; 
then} j,=— > = ~e— a") 47 (a—| 

m! (—3—4r)! (m+ 4)! 8"(4m—1) 
m=0 ” 
a 
n—7? 9~*{,F,(—4—1r, 3, n4/8)—U} dn, (52) 


0 


+ x! is written for ['(2+1), as before. 
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me (m 7 ly)! —1)! 4m 
and h, S , . = ( ; s) | a 
—, m! (—{—4r)!(m— >)! 8"(2m—1) 
7 
1—27? n{,F.(—?—}, 2, n1/8)—U dn. (53) 
0 
| The series for g. terminates when r 4m-+-2, and that for h, terminates 
when r = 4m-+1, m being a positive integer or zero. 


As regards asymptotic expansions, in addition to the solution 7? the 


equation (44) with G, put equal to zero has two solutions whose asymptotic 


expansions for large » commence with multiples of 


respectively. 


n't3 andof 7-*+!% exp (74/8) 


When z is large and positive (6),7 


Ala, b,2)~ 


nd 


h 


r 


~~ 


( 


(b 1)! ee 
Ca 
(a—1)! 
b—a)(l—a) . (b—a)(b—a+1)(1—a)(2—a) ) - 
ai “ay 
3(3r+13)/4 ! 
ys ( ) 
4/8) ~ —F =“ exp(n*/8)n-“tf1 +...) (r 4 4m+2) 
(—3—4") 
(55) 
2(37 10 4/ 5 ! ' ; 
1/8) ~ (3 = exp(n*/8)n-" {1 +... (rr A 4m+1). 
4 4"): 
(56) 
2 sii *( i)! . 4/2 ~(r+10)$] + 4 | 9 wi 
(—#—4n)! exp(7 )n pi+...§ (7 A 4m+2) (57) 
3 — 4"): 
}(3r+18 4/ a)! . 
Tp! exp(n*/8)n-“" +f] +...1 (7 A 4m+1). (58) 
4 4 . 


Exponentially large terms must not occur in the solution for f,, so when 
4m- 


Butt 


1)! 


~~ (a 


( 


| or 4m+2, g, and h, must occur in the combination 


*)! (—3—4r)!9,+-2-?(— 2)! (—i—})!h,. (59) 


b)! F(a, b,x) +(a—b)! (b—2)! 21, F(a+1—6, 2—b, 2) 


Ll)! (a b)!a aly - om a 


a(a+1—-b)  a(a+1)(a+1—b)(a+2—b) | 


“Ta 2a i 


The formula is on p. 258 of ref. (6). 
t Barnes (5), op. cit., p. 259. 
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hence 
(—%)! (—$—2)! ,A(— 4-4. 4, n*/8)+ 
8!(— 3)!(—j—}r)! 9 1 ,A(—i—-}, i, 04/8) 
J1Ipr\(Ztpr 
~~ 3-r 24(_ 31)! ( z hr)! yf of, __ (2 dtd 
2 27! 
_ 2+rn2 -r)(3+r)(1— r) (2 | r)(2—r)(6—r)(3+r)(] —r)(5—r) - 
2.4.78 2.4.6.7 a 
(r 4 4m+1 or 4m+2) (61) 

and 


+ 8(—8)! (—J— I)! yf R(— 8-4, 8 nt/8)—1idy 


7__1p)! { nr+3 —_ (2+7r)(8+71)7 a 


27 ry! 4 4 “|r 17 2(r—3) 
_(2+r)(2—r)(3+r)(lL—r)y" 
: 2.4.(r—7) = 


(2+ r)(2—r)(6— r)(3 Lr)(] —r)(5— r)nr-9 ; 


2.4.6.(r—11) —* 
+-constant . 7” (r = 4m). (62) 
When r = 4m+3 the term with a zero denominator must be replaced by 


(— 1)mg—m—5/4( m- yt . m—%)! x 


5.9.13...(4m+5).3.5.7...(2m+3) , 
A SOON tage, 


(m-+-1)! 
(—1)m+1, 2 art. (—5)! i 
which reduces to y) macead . (—4) yPlog 7. (63) 


(m+1)! 
(62) may be verified and the constant multiplier of y? found by considering 


] r r nett 4 
—— | (—s—1)!(—3—s)! (s—$—}r)! s (64) 
272 | . - 4 89(4.8 = ) 

taken round a contour consisting of the straight line from —N—}—%! 
to —-N—4+00i and the part of a circle of infinite radius to the right 0! 
the line. When r 4 4m-+-3 the constant multiple of 7? in (62) is found 


to be 2-748)! (—38—I pr)! 2, (65) | 
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When r = 4m-+-3 the term with a zero denominator in (62) must be 
omitted, and (65) replaced by 


(—1)™+1772(—3)! , We ee ae ; 
ae *< 4? | log— —§(—})—B(—H)+3(m+1)], (66) 
244(m-+-1) 8 
~~ C __ 
where (2) = - -log z!. (67) 
dz 


The multiplier of 7*log agrees with that in (63). If y is Euler’s constant 


equal to 0-5772...) 


dd ooo) 


(mm 1) l += + + sr 
2 3 m+1 
, (68) 
ay | 4) 577 4 3 log 2 y 
w(—3) 1— 2 log 2—y 


so (66) is equal to 


] )@+1yz ( 2)! 2| 4 loo > | | l 4 a (69) 
97/4 ' "| S ‘i -~ st om 6 ote : 1 “ - y 2 e ° . 
2“4(m-+- 1) 2° 3 m- 


In particular, since 





! Z 7 ' Aor o\1 4.23.7 11\! 64 1! 
i) = | ( 5) ( 4): FS 1 ( 7)! oa . 
(70) 
follows from (62), (65), and (69) that 
, 3.2. . a 15 1, 15.1.3 1 
3 2 G9 : 1) — one . a. — 
10.(4!)8 160(4!)?| 2! »? 2.3! 76 
ib. t,3.e8.3 1 Sar ” is 
— = | sore 114 log n+2 log 2+y—47—5] (71) 
3.4 7) | bs ye 
and 
— 6) 7 ? = @ 
. Fe os 5 1T }7° > Se ee : 
\—— na ee t eae i : ) a 
64(4!) 48.2¢.4!|5 ] oe n 
Aen eS 8, 958.4.0.8,0 4 
(eh. n°? 11.4! 7? 
3.1.1.3.5.7.3.1.5.9 1 ) nF we (72) 
= ep — Pe d 
15.5! ;" | 640(4!)4 
We now return to the equation for fs, which is 
Ig $n°f3 nfs 6nfs —10a, Xe 1? — fait 7°. (73) 


The general integral with a double zero at the origin is 
Is X3 * + 4ay X%(9—Y3) san(1-+}%4 hs), (74) 


Where a; is a constant. In order that exponentially large terms should 
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not occur in fs, g,; and h, must appear in (74) in the same combination 
as in (71). Hence we must have 


24 ai 2 — 
fs = 53 xT (75) 
and 
r 3. 2t al 
fs = a3 7°— ait +}$n*—h;— oye (720) 


3.2? ni = __ 9)! 1! ,4m+1 
_ > (m 7)! 4:7 | (76) 
10(4!)8 me! (—§)! (mf)! 8"(4n -1){ 





f m  3f 4 151 ,15.1.31 15.1.3.3.7 1 | 
g™~ — xy {7 — - ———— oe 
. 60G!)2 2|" 272 2.3! 78 3.4! pot T 
+ a “i [4 log y+ 2 log 2+-y—3a—5] + 
i! 
3.2%. ai 
z oH 1 Tod! 3 "| (77) 
Oo} 
Hence, from (42), a, = ———__, al. (78) 


We now turn to the equation for f,, which on substituting in (44) we 
find to be 


fi in3 fi +4n?fy —Tnf, = P,—6(a§+ 2a, 3)? — 2a¥ Xo 7 


a | % (16m?— 20m-+-3)(m—3$)! (—4)! 4” 
— way - & 
ali — m!(—§ 5)! (m—4)! 8"™(2m—1) 


a > (16m?— 12m—1)(m—$)! 4! 4”) (79) 


m! (—{)! (m+-})!8"(4m—1) J 
The general solution of (79) with a double zero at the origin is (with 


denoting a constant) 


105 


; as n? a 77! 
Ja = %q°+ a(n r) +-2(08-+ 2a ¥3)(9—9J4) — Bat “(1+ c —h}- 
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where 
ba Sm) (4)! (2m —3) gi? 
H+, m! (—3)!(m— 4)! 8"(4m-+-2)(4m-+- 3) 


.t.> (m—{)! 4! (4m—5)n™4 (81) 
40($!)® Ze am! (—$)! (m+ 3)! 8"(4m+3)(4m-+4)" 
Now ” 
d2L ‘ d { ] 


7 ’ Fy 
dy? ' dy \294 ' 


{.Bt/8)—1+4n4]- 
3.2.2? 1 ay: 

— —__. —_|, Ff, (— 8, § »*/8)—1}}, 82 
40(4!)8 alt i 4°¢ i) ) ] ( ) 

and since, from (61) and (70) 
] : = ae 3.2¢.7? | 

u —‘ , 7 Ss —_ 

ghi(—2, 4 0'/8) 40(4!)? 7 


1] $ 3.7.1.3 3.7.11.1.3.5, ) 


3] 
OTT | 


ae, 
9 > ~ s T 7 ‘ oT 4 ~ 6 ’ 
B2(E2(15 yt Was Bly" 4t yi 5! 20 j 


~ (83) 
the asymptotic expansion of L may be deduced apart from an additive 
constant and an additive multiple of 7. The resulting expression may be 
checked, and the constant and the multiple of 7 determined, by considering 
i . 4s—5)n*+4 

_ 1)! (—s—8)! (s—2)! { - ] ds (84) 
27 8*(4s-+ 3)(4s+ 4) 
taken round a contour consisting of the straight line from —N—}—ooi 
to —N—4+coi and the part of a circle of infinite radius to the right of 
the line. After some calculation it is found thatt 


7 9} 3 3 2 
a”. OTT , 
es i ~  < nf4log n+2 log 2+y—497—3!— 
9s 5 7. es © £ PY; 9 TN 57) > as Y 2 “3 
168 32(4!)8 2 32(}!)? 
3. int = Bax [x [ t. eet , Loe 1 
5 +3 ix as <a —- 
20(4!) 2(4!)? | ) 1.2 7° 2.3! n 3.4! 11 | 
(85) 
As before y denotes Euler’s constant, and use has been made of (70), of 
the first equation of (68) with m 1, and of the formulae 
(—t) t7—3 log 2—y, (—4) —2log 2—y. (86) 


The asymptotic expansions of all the terms in the expression for f, in (80) 
are now known. In order that exponentially large terms should be absent 
from the asymptotic expansion of f,, g, and h, must occur in (80) in the 
same combination as in (72). The terms containing g, and h, in (80) are 


8.2%. xz? 
¢ ‘ ¢ 52 decile: 7 
2( X9 te XY X3)94 +(} xy %e-+--— 11\3 xt hy. (87) 
7(¢!) 


' I am indebted to Mr. C. W. Jones for the correction of three errors in (85) as first 
written 
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Compare with (72), and use the value of a, from (75). We must have 





iy Po) 
D(~2_-1¢ sa | , 
2(a3 Pp Oy 3) 20462 (88) 
3 
7 
_(35—8. 93)q3 RC 
and Xs 4000) (35—8. 2*)a¥. (89) 


Substituting for the coefficients of g, and h, in (80) from (88) and (75), we 
find that 


— n* sos 
Is er Ce r) — Hat L-+ 


: .2t 7 { 77! 7.2%. 73 
4 
+ sar tlhe) agp eo} 0 


The series expansions of L, g,, and h, are given by (81), (52), and (53), so 


that the series expansion of f, is easily written down. Its asymptotic 
expansion is found from (72) and (85) to be 


, 4 P\, 2@.nt , a, [FB = 
I~ & ~ 630)" 15(g)3 27 +(2- w+ 


21.2%. 24 4a 
+(%— Anos 17 at) ton 
+! 


x4 » log y+ 
400(4! 08 1" 


1201 


ot F V7 
Tae xi "| log 2++-y—30— + sal 
77 {7 3.7 a Ab Bonsai 1 |] 
eet we arch elndh ies 
0g)A le 1” i! 7.3!» 
PASSS4181, 1 
11.4! ne J 
Sag) 3 SOS t TARTS ly 
( ')2 15 1.2! »? 2.3! 7» 3.4! oo | 
It follows from (42) that 
oo > 
a, = of| = — a | <a (92) 
9 600(4!)4| 360 





The analytical discussion of f, is now complete. The equations for f;. 
fe, etc., have not been completely solved, but some discussion of these 
equations will now be given. 


In addition to the complementary function «,?, the equation for f, has 
one complementary function whose asymptotic expansion for large positive 
values of 7 commences with a multiple of »-*+%exp(7*/8), and another 
complementary function whose asymptotic expansion commences with a 
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multiple of »’** followed by a multiple of y’-!. For the solution to be 
successful, exponentially large terms must not occur in the asymptotic 
expansion of f,, which must begin with a multiple of 7’+*. For r < 4 this 
condition is satisfied, and, moreover, the next term in the asymptotic 
expansion is a multiple of 7’*!. Assume for the moment that these state- 
ments are correct for r < n—1. Consider the asymptotic expansion of 
@,, the left-hand side of the equation for f,, in (44). The terms of highest 
degree occurring could be multiples of "+4, but it is not difficult to prove 
that the terms in »”** always cancel, and that in fact the asymptotic 
expansion of G,, begins with a multiple of n"*+*. It follows that the equation 
for f, has a particular integral, which we denote by J, whose asymptotic 
expansion begins with a multiple of »”"*?. Any particular integral may 
be expressed as the sum of J and multiples of the complementary func- 
tions. The particular integral with a double zero at the origin is indetermi- 
nate only to the extent of an additive multiple of »?, and what is required 
is that it should not involve the complementary function that is expo- 
nentially large at infinity. If the presence of the exponentially large 
complementary function can be avoided, the asymptotic expansion of the 
solution with a double zero at the origin will begin with a multiple of 
7"*% (from the other complementary function), followed by a multiple 
of 7"*! (from J). Hence by induction it will be true generally that the 
asymptotic expansion of f, will begin with a multiple of 7”**, followed by 
a multiple of "+1. 

Now /f,_, contains a term «a, _, 7”, where a,_, is a constant which is 
undetermined at that stage of the solution at which we solve for f,,_;. 


The only terms in G,, containing «,,_, arise from the terms in the expression 


for G,, which contain f,,_, or its derivatives, and it is easy to see that the 
sum of the terms in G,, containing a,_, is —(2n+4)a,a,_,7?. The corre- 
sponding term in the solution for f 


Jit 


with a double zero at the origin is 
4a, x, -,(9—g,,). Unless n = 4m-+-2, where m is a positive integer or zero, 
the asymptotic expansion of g,, involves exponentially large terms. Other 
exponentially large terms occurring in the solution for f,, can arise only 
from multiples of g, or (when n 4 4m+1) from multiples of h,. Since 
a suitable combination of /:,, and g, has an asymptotic expansion devoid 
of exponentially large terms, and this combination involves h,, unless 
n = 4m--2, it follows that when n + 4m--2 the presence of exponentially 
large terms in the solution for f 


Jn 


can always be avoided by a suitable 


choice of «,_,. In this way, when n ~ 4m+2, «,_, is determined. 
On the other hand, when n 4m+2 the series for g, in (52) 
terminates, the term of highest degree in g, being a multiple of "+. 


In such cases it is not possible to arrange for the absence of exponentially 
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large terms in f,, by a suitable choice of «,_,, and some other condition 
must be satisfied. 

At each stage, in solving for f,, the value of a,,, is fixed by (42) 
in terms of such of the P, as have occurred in the equations up to and 
including the equation for f, and of such of the a,, for r << n—1, as 
have not been determined by the conditions for the absence of exponen- 
tially large terms. 

In order to proceed farther, we consider in some detail the equations 
for f, and f,. The equation for f, may be written 
7 


fs —3n°fs +3n7fs—8nfs = —140, x4 7?—§a,P, (842, 


Jafri), (98) 
where H,(n) is a function of 7 independent of a,, a4, and P,, and the solution 
with a double zero at the origin is 


2 oF, 61 5 ¢ 
fs = 5 9? +404 x4(9—G5) ———* 8 + af k(n), (94) 
4 


where k; (7) is a function of » independent of «,, a4, «;, and P,. The presence 
of exponentially large terms in the asymptotic expansion of f; can be 
avoided by a suitable choice of a,, of the form 
_ aq = (constant )af. (95) 
The equation for f, is then found to be of the form 
mm” 136" | poe oe 
Se —37°f6 +5n‘fe—Infe P 
7 ; 2_26.2p 168, pi(,3, 7 6 9F 
= — 16a, «5 9° —f504 P, n° — fa, F, (» +7) + af Hg(7), (96) 
where H,(7) is independent of a,, «4, a;, and P,, and the solution with a 
double zero at the origin is 
2.2 Pp 
13aj P, 
11340 


gl 
45 


9 


la 


te = %g 1+ 40, x5(7 —Jg) — 7+ x8 ke(7), (97) 


where k,(7) is independent of o,, a;, «,, and P,. From (52) 
5 9 
7 
$= 1+5-. (98) 
15 1260 


Hence 


2 
2461 Ska(n). (99) 


7? 47° = 13 x? P, “ > 
11340 45 


P 9 

» = Ag 7+, Oe — 
to 67) 1 (50 15 
In order that the asymptotic expansion of f, should contain no exponen- 
tially large terms it is necessary that the asymptotic expansion of k, should 


contain no exponentially large terms, i.e. that the particular integral of 


fa —Anfe + 5n%fe—Infy = Hel) (100) 


with a double zero at the origin should contain no exponentially large 
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terms. This particular integral is indeterminate only to the extent of an 
additive multiple of »?, and since the expansion of H, in ascending powers 
of » begins with a multiple of »*, we may consider the condition to refer 
to the particular integral of which the expansion in ascending powers of 7 
begins with a multiple of n°. The condition may be given a more definite 


— > 
form. Put l=, on Y= % (101) 


Then (100) becomes 

n?z4+(6n—475)z,+ (6+3n')z, = A,(7). (102) 
A complementary function of (102) (corresponding with the complemen- 
tary function g,(7) of the equation for f,) is 


1 7 6 
Ue = a: a (103) 
: ae 5 180 
If we put Mg = Mgte and vv, = Wp, (104) 
(102) becomes 
oe . 
: [ uz n®w,exp(—7*/8)| = H, ug nt exp(— 74/8). (105) 
( 1) 


The expansion of the right-hand side of (105) in ascending powers 
of » begins with a multiple of »*. As explained, the expansion of the 
required particular integral for f, begins with a multiple of °, and hence 
the expansion of u2 y°w,exp(—7*/8) begins with a multiple of 7°. The 
required solution of (105) is therefore 
7 
uz 1®w, exp(— 74/8) H, Ug, ntexp(— 74/8) dn. (106) 
0 
If the asymptotic expansion of f, contains no exponentially large terms, 
the asymptotic expansion of w, contains no exponentially large terms, 
and the left-hand side of (106) +0 when »-o. Hence the required 
condition is 


x 


( H, Ue ntexp( -»!4 8) dn 0, (107) 
0 
: . 76 70 
Le. H,\ 7? ——+ a0 exp(—7*/8)dyn = 0. (108) 


0 
The expression for H, is easily found from the expression for G, in (46), but 
[ have not been able to evaluate analytically the integral on the right 
of (108). The matter has been further considered by Mr. C. W. Jones, who 
finds numerically that the integral condition is satisfied to the accuracy of 

his computations. 
If we assume for the present that the condition (108) is satisfied, we 
may proceed to consider the equations for f,, fs, etc., in the same way as 
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we considered the equations for f; and f,. From the condition for the 
absence of exponentially large terms in the asymptotic expansion of f,, 
v, is determined in terms of a,, «;, and P,. Similarly «, is determined from 
the equation for f, in terms of «,, a;, and P,,f and so on until we come to 
the equation for f,). In the right-hand side, Gio, of the equation for f,,, 
vy occurs only in the term — 24a, a, ?, and the corresponding part of the 
particular integral with a double zero at the origin is 4a, a9(7—g,9). Since 
the series for g,) terminates, the term of highest degree being a multiple 
of 7’, this part of the particular integral contains no exponentially large 
terms, and a, will not be determined from the condition that the asymp- 
totic expansion of f,) should contain no exponentially large terms. Gy also 
contains a multiple of a2? and certain very complicated terms linear in 
x5, in particular a term in afa;. (It further contains in particular compli- 
cated terms that are multiples of aj° and af P,.) The part of the particular 
integral that corresponds with the multiple of a2? in Gi») will be a 
multiple of a2(7—g,9), and will not be exponentially large at infinity. 
Unless, therefore, the part of the particular integral corresponding with 
those complicated terms in G,,. which are linear in «a, also fails to become 
exponentially large at infinity, «, will be determined from the condition 
that f,) as a whole should not be exponentially large, which condition it 
will thus always be possible to fulfil. In other words, it seems probable 
(though I have not proved it) that «,; is determined from the equation for 
fio, % from the equation for f,,, and so on. If this is correct, then only 
x, remains undetermined among the «’s, and therefore only a, among 
the a’s. As explained in the introduction, a, is then probably determined 
by the condition u,—>1 as y, > at 2, = 0, in which case the whole 
solution is determined at separation, and probably what we have con- 
sidered is an asymptotic solution at and near separation, applying (in the 
non-dimensional form here considered) to all cases in which the P’s are 
the same. 


3. The solution upstream for large values of y,/x} 

Our first aim in this section is to exhibit the form of the solution for 
large values of y,/x}. By carrying through the suggested calculation in 
some detail we also provide a check on part of the work in the preceding 
section. 

- 


; P , ; P? 
+ The right-hand side, Gg, of the equation for f, contains P, and a term, — = (at+t), 


in P?. The corresponding parts of the particular integral with a double zero at the origin 


are, however, : 7 
A(t-z ia) oe ue 
*"\6 315° 31185 1\630 ' 155925 


and contain no exponentially large terms, so that the expression for «, will involve neither 
P, nor P}. 
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he According to the solution in the preceding section, is found in the 
I ~ 1 
eo form of the series (24), where 
ym 73 x2 . 
’ / 9 » [- 
t Jo ;? F oh fo = 7-7, (109) 
oO 6 15 
10> y. is given by (75), and for large positive values of 7 
he : ' “a ‘ , , 
1ce fs ~ As °+C  n*+ Es 4? log n+ Es 4?+ Fs + G3+Kgn?+..., (110) 
ple where 
—_ : , . 7 ) 
ge é i 3 ? 11\2 vy ( 3 3 vi E, 11\2 “e 
O0(4 ” (3!)° | 
1p : 
7 iT ‘ 77 - 
lso E. 7 a3! 2 log 2 to —5-+- q(35—8 2)\, { (111) 
| 4(4!)2 “| 100(4!) | | 
i 
4.2% a? . , 7 
: . ‘ 3 y g 3 i 3 | 
pli F, Tl eee Gs 3}; K, = r m3) 
lat 4) 3(q') J 
and 
a 
- ” nl 3 ’ 2 y/ , - 1 
oe ly ™ A,? tOU,y dD, n E47 1 Fy» : G7 log n +G,n+Kyn 14. ..., 
ith Ges 
where 
me { a ; P. T . ) 
i A X ‘ 9 = — on ew Le » 
a é‘ 63 150(4!)4] 7 630 . 5(4!)2 7 
1 it 5 ei 
. wa > » 2 
; D = 7 Af E : 2 16 in 4 
le 4 15 11\3 1> “4 , 3 11\4 Xy> 
f 5(4') 6 : 10(4!) 
or | 
' 91. 292. x4 
nly | F, : “= at (a has not been calculated), (113) 
Pt 400(4!) 
oa , dcr T 77? 
e | Ai 4 G 419 a .2 = 
- = Gp — Gaqyaat| 2los2+y—-de-34 SI, 
ole } . 7 4° 
7 a 4 
on K, = oped 
the o"4:) J 
are Now substitute the expressions for fo, f,, f. and the asymptotic formulae 
lor f; and f, into the series (24) for %,, substitute for » from 2*én = y,, 
and rearrange in powers of € and logé. The work is purely formal, and 
P no justification is attempted. The result, which can hold only for large 
‘ ) 
¥ values of Y,/v7, 18 
1 in usu 
5 Ma lg 2 Z - 
ling ub FS an ae Oe SyS 4 —Sy7 st 
1 6 30/1 Tks 4 4 
1), &2(2 M1 Yi +2 C5 yi +3C,y2+...)4 
15 “ee 1 > 1 | 
a + 8( hay 92+ 2-1 Dy yt +...) + 
+ &4{[24H,— 2-1 Es log 2 y+ 2! E y? log y, + Ey y+...}+ 
E4 log €(—2* Hy y?+...)+ 
ither 


5(2F; y, + 257, y2+... 


Sc 





(114) 
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We are therefore led to assume, as a form valid for sufficiently large values 
of y,/24, 


By = XolYW)+E?xol(Yi) +E x9(Y1) + E4xa(Yr) + (4 log E)X4(y) +E x5(¥i1)+--, 


and the initial terms of the expansions of the x in powers of y, are found 
by comparing with (114). In fact, since the form (115) is not valid unless 
y,/x¢ is sufficiently large, the boundary conditions at y, = 0 cannot be 
applied, and the solution found must be joined to the solution in the 
preceding section by using the first few terms of the expansions of the y’s 
as given by (114). On the other hand, the solution in this section should 
satisfy the condition u, > U, as y, > ©. 
From (115) 


lu, = Oy = xo tOxet Syst ety t (Flog €)x, +E? y¥5+...; (116) 
OY, 
i ] 
y= _ Hy —— {2x X2° +- 3Ex5 i i E7(4y4+X4)- r. 4(£" lo fE)X, +5& Xs woofs (117) 
c “a 4e° 


We substitute in the reduced equation of motion (16) (with é@p,/éx, as in 
(17)), multiply by 4€?, equate coefficients of €°, €, €*, €@log €, and é°, and 
obtain the following equations: 
Xo X2—Xo X2 = 9, 
Xo Xs— Xo Xa = 9, 
Xo(4x4a+Xa)—Xo(4xa+ Xa) = 4+-4x9 +2(x2?— xe x2); (120 
X0 Xa—Xo Xa = 9, ( 
5(xX0Xs—X0X5) = 5X2 X3— pan “> Xs: ( 
Now Xo = dyi— tay tyi+ ay 4 4 yi +.. (123) 
and is, of course, the value of uw, at € = 0, Mio 
byit+asyita,yit+ 
as may be verified from the values of A;, Ay, a3, G4, 5, Qs. 


The solution of (118) is 
X2 = (constant)x. (124) 


The constant multiplier is determined by making the coefficient of Yi 
in the expansion of xy, equal to 2'a,. Hence 

hee 2a, Xo- (125) 

(It may now be verified that the coefficients of y} and y? in the expansion 

of x, are the same as in (114), and the values of C, and C, are thus checked.) 

Similarly from (119) and (121) it follows that x, and x, are constant 
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multiples of yj, and the multipliers are determined from the coefficients 
of y? by comparison with (114). In these multipliers the values of «a, and 
FE; are substituted from (75) and (111), and it is found that 





9% zi : 
mo Bays 1X (126) 
a: 
P 2ta , 
X4 (q1yeo Xo" (127) 
4 


(The coefficient of y} in the expansion of x, may now be verified, and the 
value of D, so checked.) 

In (120) the terms in ¥, may be omitted, since they cancel by (121), 
and the value of y, may be substituted from (125). The equation for 
y, is thus seen to be 

(28) = x08 x5 Aa xa?— x0 x8} (128) 
dy1\xo 
The expression on the right of (128) may be expanded in a series of ascending 
powers of y,; the first terms are 
- “2+ of -2104,)+.. (129) 
Hence the solution of (128) is 


eke xo — 49i(xo"— XoXo) , 120.24. As| 
J | Xo Y J 


—120.2*. A, xg log y,+(constant)xy. (130) 


dy, — 


© 
0 


The constant multiplier of yj is found from the coefficient of y}? to be 
equal to 


n 
100(4!)4 


5) 


NE, 2% log2 = —" Hy jr —5 4 (35—8.2!)|. (131) 


2(4!)* 
The coefficients of y? logy, and y? may now be verified, and the values 
of E; and E, so checked.) 
When we substitute for y, and x, from (125) and (126), the equation 
122) for x; reduces to 
8. 2t.7 


Sq 8d (XO — Xo x0) (132) 





XoX5— X0X5 = 


rate 


with the solution 
S.2*.a 


Xs —-— a} Xo +-(constant) x9. (133) 
o(4z:)° 
from the coefficient of y? the constant multiplier is seen to be 2#F,, but 


it cannot be fully determined from results already found, since the 


expression for F, in (113) contains a,, which has not been calculated. 
5092.1 


F 
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(The coefficient of y, in the expansion of x, is determined without a 
knowledge of the constant multiplier, and is easily seen to be 2F;, as in 
(114). The value of F, is thus checked.) 

It remains to consider the condition that the velocity should pass 
smoothly over into the velocity of the main stream. As explained, we 
must suppose that this condition is satisfied at x7, = 0; as far as our 
formulae go, u, is given by (116), so if the condition is satisfied at x, = 0, 
Xo > 1, and the second and higher derivatives > 0, as y,->0. From 
(125), (126), (127), and (133) it follows at once that x5, v3, %4, x; and their 
derivatives all > 0 as y, >. As regards x4, we suppose that, as y, >, 
Xo > 0 more rapidly than y;'; it may then be proved that as y, > 0, y, is 
asymptotically equal to a multiple of y,, plus a constant, plus terms which 
+> 0, so that y, > © as y, > 0; but x5 x, > 0, and we may then show from 
(118) that y,—> 1 and the second and higher derivatives > 0, as y, >. 
We thus check that, as far as (116) goes, u, > U, (since 4 = x, and 
U, = 1, dU,/dx, = 1 at x, = 0), while @u,/ey, and higher derivatives — 0. 


4. Discussion of the solution upstream 

The solution in the preceding section holds only for large values of 7. 
The solution in §2 applies for small or moderate values of 7; it is valid 
also for large values of y (the asymptotic expansions of f, and f, being 
used), but is then useful only for small values of y,, since in such a case 
it is essentially equivalent to the solution in §3 with the x’s expanded in 
series and with only a few terms in each expansion known. The results 
that can be obtained from §2 when 7 is large may therefore be more 
advantageously obtained from §3 after we have obtained the formulae 
for the coefficients a, in the expansion, in powers of y,, of the value of 
u, at x, = 0. (It should also be remembered that the asymptotic formulae 
of § 2 were necessary for the completion of the solution in § 3.) 

Thus the formulae in §§2 and 3 are useful in different regions. In 
particular, the formulae of §2 will be useful for studying the values at 
y, = 0, x, # 0 of the derivatives of u, with respect to y,, etc., whereas 
the formulae of §3 will be useful for studying the nature of the solution 

), » ~ nla, 


The formulae of § 2 show that lim (é”u, /é2 
The formulae of § 2 sl that lim (é"u,/éy7 


z1—0 


for n = 5 and 6. For any non-zero value of y,, however, @”"u,/éy{ 1s 


at x, 0, ¥, F ( —_ 


continuous at 2, = 0. On the other hand, the formulae of §3 show that 
v, (and éu,/éx,) are infinite at x, = 0. 

Because of the singularity at 2, = 0 the usual assumptions of boundary- 
layer theory are invalid at 2, = 0 and in the immediate neighbourhood. 
Nevertheless, the mathematical result that v, is infinite may be taken to 

1 . 
indicate that large cross-velocities are to be expected at separation; 
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otherwise the assumptions of boundary-layer theory would not break 
down. 


5. The solution downstream 

In considering the motion downstream from separation we reduce the 
equations of motion and continuity to non-dimensional form by the same 
substitutions as before, except that in place of x, we use 


, 


ra —2, = (x—2,)/l, (134) 
so that x, is positive downstream. The governing equations are obtained 
by replacing x, by —2, in (16) and (17), and in place of (23) and (24) we 
write be ' oir on 

E x}, n' = y,/2*2’t, (135) 
py = 283 Fy(n’)+€'F(n')+€2F(n’)+...]. (136) 
We obtain differential equations for the F in the same way as before; 
moreover, the boundary conditions are the s»me, for the conditions ys, = 0 


and u, = 0 at y, 0 lead to 


F.(0) F’(0)=90 (r =0,1,2....) (137) 
and the condition that lim w, is given by (21) leads to 
£0 
lim —~, = 2#a,,. (r = 0,1, 2....). (138) 
n> « 7) san 
The equation for F, is 
F" +3F, F7—2F? = 1, (139) 
and, since a, = }, the solution is the same as for fy, namely, 
F, = 73/6. (140) 
The equation for F, is 
FY +4F, —$n’*F\+47’F, = 0, (141) 
and, since a, = 0, the solution is 
F, = B, 9", (142) 
where 8, is a constant. The equation for F, now becomes 
Fy +47 F3—3n'2F3+57'F, = 46? y”, (143) 


and the general solution for F’, with a double zero at the origin is 
a Q /21 192,45 
F, = Bg “+-35P1'"s (144) 
where 8, is a constant. Hence, from (138), 
- i K 
a, = +63. (145) 


2 . ~ . . . . . 
But, from (50), @, is zero (in which case a, is zero) or negative. If a, 
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is negative, there is no real solution downstream of separation. Hence 
there is no real solution downstream of separation unless a, = 0. 


6. The special case a, = 0. The solution without singularities 

If we return to §2, and consider the motion upstream of separation 
when a, = 0, the equations for the f are easily integrated, and we may 
show that a, = a4 = a; = ag = ag = a = etc. = 0; the whole solution 
is free from singularities and may be verified by expanding y, in a double 
power series in x, and y,. Moreover, since the solution is free from 
singularities it will hold also downstream. It is not possible, however, by 
the methods used, to consider if there are pressure distributions for which 
U,>lasy,>0 atx, = 0. 

The most interesting special case of this solution is that in which 














YX = a =... = 0, when it reduces to one that is easily found inde- 
pendently, namely (with x, = —z,), 
5 Op,\? Pp, 9 Pr op,\" 
oe 1 dp, Mi 1 ap, Yi dx] oa;® dx, \oax'? a 
1 2 ax, 3 ' 360 ax? 7 453600 De 


(146) 


so the expansion of w, in powers of y, contains only multiples of yj". 
This is a solution in which (éu/éy),.. = 0 for all x, but in order that it 
should be valid it is necessary that @p,/éx, should be chosen so that 
u, > 1 when y, >© at x, = 0. Included in this solution is that special 
case of the solution discovered by Falkner and Skan (4) for which 
(6u/éy),-9 = 9 at all values of x. In the case considered by Falkner and 
Skan the velocity distributions at different values of x are similar; if more 
general solutions of the type shown in (146) exist, the velocity distributions 
at different values of x will not be similar in the general case. 

It is a fairly straightforward matter to check that the known solution 
for U = cx™, when m has the appropriate value, agrees with (146) as far 
as that equation goes; but the value of m is determined from the condition 
u— U as yoo and cannot be found by the methods used here. Since 
the velocity distributions at different values of x are similar, the appro- 
priate value of m may be found from the solutions of an ordinary differen- 
tial equation, and has been so found by Hartree (5).f No such method is 
available in the general case. Meanwhile the formulae of Falkner and 
Skan, when (0u/@y),-9 = 0, have been fitted as a very special case into 
the formulae of this section, so far as those formulae go. 

+ For negative values of m the solution of the equation with the conditions y = 0, u = 9 


at y = 0 and u/U — 1 as y > is not unique, but may be made unique by requiring that 
1—u/U shall be positive and shall — 0 exponentially as y — oo. 
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NOTE ON THE VELOCITY AND ‘TEMPERATURE 
DISTRIBUTIONS ATTAINED WITH SUCTION ON 
A FLAT PLATE OF INFINITE EXTENT IN 
COMPRESSIBLE FLOWT 


By A. D. YOUNG (Department of Aerodynamics, 
the College of Aeronautics, Cranfield) 


[Received 4 November 1947] 


SUMMARY 

The problem considered by Griffith and Meredith (1) for incompressible flow is here 
considered for compressible flow of a perfect gas with constant specific heats, it 
being assumed that there is no heat transfer by conduction at the plate. Essentially, 
the method consists in establishing a correspondence between the velocity and 
temperature profiles for incompressible flow and those for compressible flow, the 
lateral ordinates being scaled by factors which are functions of the ordinates and 
of Mach number. 

The results of calculations covering a range of Mach numbers up to 5-0 are shown 
in Figs. 1 and 2. 


1. Notation 
x distance measured parallel to the plate in direction of main stream 
upstream of plate 
y distance measured normal to the plate from the surface of the 
plate 
u velocity component in x-direction 


v velocity component in y-direction 

p density 

T temperature 

ph ~—s coefficient of viscosity 

k thermal conductivity 

¢, specific heat at constant volume (assumed constant) 
c, specific heat at constant pressure (assumed constant) 
o wc,/k (Prandtl number, assumed constant) 


J mechanical equivalent of heat 

c,,/C, (assumed constant) 

i Jc, T (enthalpy) 

T Pr (shear stress) 
dy 

suffix 1 refers to quantities measured at large normal distances from the 
plate (y > «) 


ie 3 


{+ This paper is substantially Report No. 8 of the College of Aeronautics. 
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suffix w refers to quantities measured at the plate 


w defined by() = (7) 
My T, 


Yr) 
ui] 

“" 

7 

L 

ft (' / dy 

J \B 

0 
y ur | (a, is the speed of sound in the main 
ta, ON \(y—Ie,T} stream) 
6 ii, 


b \y l ) Mi. 


2. Introduction 

The solution due to Griffith and Meredith (1) of the velocity distribution 
attained with suction on a flat plate of infinite extent in incompressible 
flow is of special interest, since it is a solution of the general equations of 
motion and does not depend on the usual assumptions of boundary layer 
theory. The corresponding problem for compressible flow is by no means 
as simple in its most general form. However, if the usual assumptions 
of boundary layer theory are made, it permits of an exact solution which 
is easily obtained. In the following it is assumed that the gas is perfect 
and that the specific heats and Prandtl number are constant. 


‘ 


3. Analysis 
The equation of motion in the boundary layer of a flat plate at zero 
incidence in steady compressible flow is 


tt ant < (ner). (1) 
Ox oy poy oy 


The equation of continuity is 


< (pu) + —(pv) = 0. (2) 
Cx Cy 
The energy equation is 
Jc, Uu oa + Jc, v- r == J od Ph ye a a (3) 
~ O82 oy pcy\ oy p \cy 


We are interested in the problem of the final velocity and temperature 
profiles far downstream from the leading edge of the plate when all 
quantities are independent of z. 
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Hence, the above equations become 


yuu = d du 4 
r dy dy May : (4) 
pv = const. = p, 2, (5) 
(di Audi) | (ae? ' 
p dy  dy\o dy . dy} ’ 48) 


Prandtl’s number, assumed constant). 


+ Cc 
where 7 = Jc, T and o = ee ( 


The gas equation leads to 


p, 7 re 


It will be assumed that the variation of » with T is given by 


mw _ (ty _()" ) 
Py ia (7) (*) ®) 


where w = const. For air at normal temperatures w is about 0-76, but it 
increases slightly with 7. 
The boundary conditions are 
Sa a 
uU= Uy, p=—py, Y=, t=h, 7; =0 at y=, 
y 
di 


u=0, —= 
dy 


0 a y=, 


if no heat transfer by conduction is assumed to occur at the plate. 
If in equation (6) we change the independent variable from y to 1, 


fe du . : eet ; ‘ 
writing 7(w) = a i = iu) and eliminate pv by means of equation (4), 
y 
we obtain : as 
dr di d*1 
1—o)— —+7ri/—-+o]=0 9) 
( i du (i 1 ) \ 
From (4) and (5) i d ’ 
‘TO é 0 ; —_— = i= lan 
du P Pi% 
and hence T= p,v,u+C, 


where C, is a constant. 
If we write 7,, for the value of 7 at the wall, C, = 7,,. Further, since 
- — ‘_— ’ ml TY, 4 
t = 0, when u = Wy, C, = —p,v,u,. Therefore 
T= —)p V4 (u,—%) ) 


= Pp, Vy U+T,y. 


(10) 
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Equation (9) can then be written 


a 7) Fo9+ 2) =0, (11) 


a) 
du 
with the solution 


ou 2 2 o 
a (12) 
2(2—a)| Uy o u,} | 


satisfying the conditions 7 
di 
du 


At the wall, where x 0, 


i,, When wu = u,, and 


-Q0, when u= 0. 


uz 


‘ ° L] 
es (13) 


- 


and hence the total energy at the wall differs from that in the main stream 
only by the quantity (v?,—v7j). 


: ; du 
From (10) we have, since 7 p—, 
dy 
“d 
u 1 —exp( ‘ { %). (14) 
Uy J 
: 
Let n > = tL 
Vy 
and let —d{=dyn ae 
od 
7 w ry w 
1 1 == d 1 3 15 
tala) = (7 a 
with € = 0, when 7 = 0. 
Then, from (14), > ian 1—exp(), (16) 
Uy 
: ; : our | 2 | - 
and from (12), ,—1 exp(2¢) ——exp(of)). (17) 
2(2—o) | a 
Writing 6 = i/i,, b = (y—1)M?, then 
6—] o 
= — F 
ip — (¢), 
: (18) 
9 
where F(¢) = —exp(ol)—exp(2¢) 


Co 

From (16) and (18) we can express u/u, and (@—1)/46 as functions of ¢ 
only, independent of Mach number. To derive the actual velocity and 
temperature distributions for any given Mach number we need to evaluate 
the relation between { and 7 (or y) given by (15). 
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c t 
_ i\” bo 
From (15), —7 = { (*) dl = Ip ret) dl. (19) 
0 
In general, the integral on the right-hand side of (19) must be evaluated 
either numerically or graphically, giving » as a function of ¢ and M.. 
Since v, is negative, only negative values of € need be considered and it 
will be found that values of |¢| greater than 10 may be ignored. Having 
determined 7 (or y) for a comprehensive range of values of ¢ and J, we 
can then, for each Mach number, replot w/w, and (@—1)/4b as functions 
of 7, using the basic (or incompressible) profiles given by (16) and (18).+ 
For the special case w = 1-0, (19) can be integrated outright to give 





_ b ao {2 exp(2Z) 2 1 ; 
=—_ ots 29 \g2ex Pld) — 9 —sat5 . (20) 


4. Calculations and results 
The velocity and temperature distributions have been calculated for 
w = 0-76 and M, = 0, 1-0, 2-0, 3-0, 4-0, and 5-0, o being taken as 0-72. 
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For comparison, calculations have also been made for w = 1-0 and 
M, = 1-0, 3-0, and 5-0. The resulting velocity distributions as functions 
of 7 are shown in Fig. 1, and the corresponding temperature distributions 


t This process of establishing a transformation of the lateral ordinate y, which converts 
the temperature and velocity profiles for incompressible flow to those for compressible flow, 
was used by Hantsche and Wendt in ref. 2. They then applied it to the boundary layer on 
a flat plate in compressible flow without suction for the special case where w = 1-0. How- 
ever, it seems capable of much wider application, and it is hoped to use it for more general 
problems of the boundary layer on a finite flat plate both with and without suction in 
compressible flow. 
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75 
are shown in Fig. 2. It will be noted that there is a thickening of the 
velocity and temperature boundary layer with increase of Mach number, 
ind this process is enhanced by an increase of w. 
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SUPERSONIC FLOW PAST SLENDER POINTED BODIES 
OF REVOLUTION AT YAW 


By M. J. LIGHTHILL (Department of Mathematics, 
The University, Manchester) 


[Received 10 November 1947] 


SUMMARY 


The effect of a slender body of revolution of thickness ¢ and arbitrary shape, held 
with its axis at an angle % to a uniform supersonic stream, is considered. By use 
of the exact equations of isentropic motion an expression for the pressure coefficient 
is obtained accurate enough to give the lift and moment coefficients as sums of 
terms of orders uf, 3, Zt?, and Wt@logt-1, and the drag coefficient as a sum of terms 
of orders ¢?, #logt—1, and y?, terms of smaller order (only) being neglected in each 
case. 


1. Introduction 


APPROXIMATIONS to the lift and moment on a slender pointed body of 
revolution whose axis is inclined at a small angle y to a uniform supersonic 
stream were obtained by Tsien (1), by solving the ‘linearized equation 
of motion’ both approximately and accurately. For a shell his first 
approximation can be written C, = 2%, and this is correct if terms of 
order ¥° and yd? (¢ the maximum thickness) are neglected. In this paper, 
by use of the exact equations of isentropic flow, additional terms in (C, 
of orders 7° and yt? are obtained, so that the new expression is correct if 
terms in 7°, £t?, and wt* are neglected. The effect of the change of entropy 
through the shock-wave will only be felt in these latter orders of small 
quantities (this remark will be supported below, § 5): hence it is legitimate 
to use the equations of isentropic flow for our purpose. In addition to the 
lift, the pressure distribution on the body and the moment of the forces 
about the nose are found, as well as the extra drag due to yaw. 

It is found that to obtain better approximations than C;,, = 2x the use 
of the linearized equation is inadequate. In a worked example the varia- 
tion of the lift coefficient with the thickness of the body and with the 
Mach number is compared with that given by the exact solution of the 
linearized equation. 

The principal application will be to shells and other supersonic pro- 
jectiles. In this paper the meridian section of the body is assumed to have 
continuous slope, in fact to be an analytic arc; a partial extension to 
curves made up of a number of analytic ares is given in a later paper (2). 
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SUPERSONIC FLOW PAST SLENDER POINTED BODIES 
2. Mathematical formulation of the problem 
In cylindrical polar coordinates (r,@,x) a uniform stream making an 
angle %& with the axis has velocity potential Uxcos~+U sinw#.rcos 6. 
Here @ is measured from the line which is the projection of the velocity 
vector on to a plane perpendicular to the axis. If a body of revolution is 
held fixed on the axis, the resulting flow, which as an approximation when 
the body is slender we assume to be irrotational (in the region ahead of 
the flat base if the body has one), will have a velocity potential ¢ which 
we will assume expansible in powers of & as follows: 
¢ = Uxcosy%+ U sind.r cos 0+-do(2, r) +h, (x, r)cos 6+ 
+-f,(x,r, 0)+d,(7,7,6)+.... (1) 
It will be found that with a potential of this form the equation of motion 
and boundary conditions can be satisfied. 
Now we have for the square of the velocity (using the suffix-derivative 
notation) 
@ = d2+¢2+r-%43 
(U cos b+$ort+$rz 008 0-+9%bo,+ Yat ---)*+ 
-(U sings cos 0+ dy,+d,, cos 0-+- bdo, +-fFh5,+...)?+ 
-r-2(— U sings.r sin 0@—pd, sin 0+-*ho9+...)? 
U?2+Q5+uQ, cos 6-4 b?V.+¥7Q3+..., (2) 
where Q, = 2U¢,,+43,+¢3,. ) 
Q, = 2(U + $o2)b12+2(U +417) For: 
QV» l Pox +2(U-4 Pox hor T $i, cos? +- 2dor pot (3) 
+ (2U¢4,,+¢4?,)cos?0+ (2Ur-14, +r-*g?)sin?6, 
Qs; 2(l I+ pox) P31 T 2¢y, COS (¢2,—$U)+ 2bo(b3-— gU cos 8)+ 
| 2(U +¢4,,)bo, cos O— 2(Ur+-¢4,)bogr-*sin 6. J 
The equation of motion is 2a*V*¢4 = (Vq?)(V¢), where 
+4(y—1 )(U?—¢q?); 
this can be written (if %* be neglected) in the form 
2{a2 — L(y ] (Yo . JQ), cos 6-4 b*Q, - WF Qs)} x 
< {Vb t+ Y(V7d, —1*d; cos 8+ PPV ha + PPV "hs)} 
ort uh, COS 6 W*bo,+ds,) x 
X (Qo, +¥Qj, cos O+-Y?Q2,+YQs,) + 
-(U cos&+do,+ud,, cos 0 +-*bo,+*¢3,) x 
x (Qor+hQ1, 008 O+- Wi? Qo, + Oa2) + 
Lr—2(— U sins .r sin 0—w d, sin 0+-*do9)(—bQ, sin O+Y?Qo9). (4) 





a” ao 


- (Usind’cosé+¢d4 
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Equating coefficients we deduce that 
2{ag— (y—-1) Qo} V"bo = Por Vor +(U + box) Gor ( 
2{ag—3(y—1)Qg}(V? 4— r*$1)—(y—1)Q:V*do 
(O+$b1r) Qt bor V1rt(O +42) Qi12+br2 Qozs (6) 
2{as—3(y—1) _ *b2—(y—1)Q,(V*d, —1-*4, )cos?8 —(y—1)Q,V7bq 
= bor Vor + (U +417) Q1, 0080+ fo Qo, +(U +2) Qort+ 12 Yi, 0080+ 
+(Por—3U)Qor+r-*(Ur+¢,)Q, sin’, (7) 
2{a5—4(y—1) Po} V*b3—(y—1)Q, cos 8 V°¢.—(y—1)Q, cos 6(V*4, —1-*,) — 
—(y—1)Q3V*do 
Por Yar + (U +447) 2, COS O+- ho, Q;, C08 8+ ($3,—4U cos 8)Qo,+ 
-(O +d or) Qs2+$12 Yor C08 8+ ($2,—3U)Q1, 608 O-+-b3, Qor— 
—Q, dog7*sin 6—(Ur+-¢,)r-*Qagsin 8. (8) 


The boundary condition is that ¢,/¢,. = R'(x) when r = R(x) if this latter 
is the equation of the meridian section. The condition becomes 


U sin cos 6+-¢5,+- bd, cos 0+-Y?d,,+Y°dbs,+... 
3 R'(x)\ U cos b+ Port bby, COs 0+ do, + Ps, + ves)s 
Hence the boundary conditions appropriate to equations (5), (6), (7), and 
(8) are . 


or 


dor = R'(x)(U+9o,), (54a) 

U+¢y, sonal (6a) 

por '(x)(—3U + $2,), (7a) 

—}U cos0+¢45, = R’ poly (8a) 

all on r = R(x). There is also a condition, which we call ‘C’, that, at 


infinity, the potentials 4; shall represent disturbances travelling in such 
a way that x increases with r. 

When the ¢; have been determined the pressure is obtained by the 
relation 


y a2\vly-D pmb 
« = ‘| =. l sie 3 > (U —q"*) 
Po i) 2019 


=1455(0—-9) +2 (0—-@ P+... (9) 








2az 8aé 
or 
1 _. P—Po ones 1 U? Nests. “ 
( ———- == —_— +4)? om 
- 4p U' ' U? ve ( - ) bi 
Q LQ cos6 w¥?Q, WFQ, , ae ae 9 
= ! <a a Ta t AMPU-*(Qo+ YQ, cos 8+ YQ. + PQs)", 


(10) 
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neglecting 4 and (1—gq?/U?); and this can be written 

C, = (—U?Q,+4M?U *Q§) +-(— U-Q, + 4£M?U4Q, Q, cos 0+ 
t-?(— U-2Q, + 4 M?U-4Q, Q.+4M?U-*Q? cos?) + 
1 3(— U-2Q, + 3 M2U-4Q, Qs +}. M2U-4Q, Q,cos 8). (11) 


The force on the body along the axis is 


27 l 
D, = 4p)U? | dé | R(x)R'(x)C,, dx (12) 
0 0 
ifthe body extends from 2 = 0 to x = I: to this a base-drag D, must be 


idded when there is a flat base. The force perpendicular to the axis 
which by symmetry is in the direction of the stream-component) is 


Qa l 
L, 4p, U? ( cos 6 dé | R(x)C, dx. (13) 
0 0 
The true lift (perpendicular to the stream) is then 
L L, cos s—(D,+ D,)sin x, (14) 
and the true drag (parallel to the stream) is 
D (D, + D,)cos b+- L, sin fs. (15) 


The moment of the aerodynamic forces about an axis through the nose 
perpendicular both to the stream-direction and to the axis of the body 
is (in the sense tending to decrease the angle of yaw) 
2a ! an F 
M -$py U*| | cos6 dé | xR(x)C, da cos 6 dé | RY(e)R'(e)C, de) 
‘0 0 0 0 (16) 
the two terms being the moments due to forces perpendicular and parallel 
to the axis respectively. The expressions (12) to (16) are all exact. 


3. First approximate solution for ¢o, ¢,, 40,5 

We now solve the equations near the body to a first approximation. 
Ift is the maximum thickness of the body, we write ‘O(t™)’ henceforth 
to mean ‘ O(t™log*t-1), for some k, when r O(t)’. 

Then equation (5) can be written V7¢, = O(t*), to which a solution is 
easily verified to be ¢, = K(x)+A(x)logr+O(t*), where K and A are 
O(#). The boundary condition (5a) becomes A R-4 U R’ + O(#), so that 
A RR’. To find K we observe that the linearized form of equation (5), 
that is, V274, = M*¢,,.., must have a solution beginning 
do K+A logr+ O(r*log r-). 


In Heaviside notation, with p- | dx and « = ,(M*—1), this linearized 


equation becomes . 


2 apy = 0, (17) 


d*4, | dd 
dr2 'r dr 
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since dy = O¢,)/ex = 0 at x = 0. By condition ‘C’ and the asymptotic 
forms of the Bessel functions Ky and J, the solution of (17) that we require 
must have the = Ky(apr)f(x); and for small r this is approximately 
{log(2/apr)—y}f(x),f or (using the Heaviside representation of a Faltung 
integral and the fact that logx = —log p—y) 

do = [ log -—y) df(y). (18) 

F ef 

This is K+A logr if and only if A = —f(x). Hence K is given in terms 
of A by the relation - 


K(x) = | oe AA ty) dy- (19) 
0 


Having obtained the relation ¢, = K+ A logr+ O(t*), we see from (3) that 
Qo = 2U(K'+A' logr) + a+ O(t*). (20) 


Similarly equation (6) can be written V2¢,—r-*¢, = O(t), to which a 
solution is ¢, = Br-1+ Er+ O(t8). The boundary condition (6a) is satisfied 
if O(t?) be neglected - taking B = UR?*, E = O(t?). E will be found in 
§ 4 by a process similar to that by which K was found: at present we need 
merely observe that ¢, = Br-1+ O(#') and so by (3), 

Q, =2u2 +2(U- 2)" +00) i ecto 
r rr} r r 





+ 0(#), (21) 


since 2A = B’. We can also write, if we assume that ¢, = O(#?), 


OTT 2 2 
Q. = ( _ 20 5 -{. vi) 0080+ (2 U - “> a] sin?6-+ O(i?) 
r pel 


r2 


= nahn Pm 20+ O(#?). (22) 


rt r2 


Equation (7) can now be written 


2a2V%4, — 1 = B 0520) +-(v—B)(—SUA 4 S42), eile: 
r4 





, r? r2 
LU Gece cos 20) +9-+(Ur +7) (A sin’#+0(€) 
r4 ia r r r 
amit tos® + cap tios® _ 100 = +0 #2), (23) 
yD r 


+ Here y represents Euler’s constant. This will not cause confusion since in the remainder 
of the paper the adiabatic index will not occur. 
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A solution is 


1 la 9A q 9 ans 9 - 

.4+1 cos 26 : 3—2 cos 20 log? 

2a? do AB: UAB = ae 
- rm) 


5U2A cos 26+ Hr-cos 20+ Llogr+N+O(t*), (24) 


where H, L, N are arbitrary functions of x. The boundary condition 


7a) becomes 
4+4cos260 .,, 3—2cos26log R+cos20 2Hcos26, L 
eo R Re TR 
- —az U R’'+-O(#*), (25) 


from which we deduce 


H = 4AB?R-+ UAB(2log R—1) 
nd L a2 URR’4+6UABR-2—2A B2R-4, 


The value of NV could be found in a similar manner to that in which K was 
found, but this is not required. We can rewrite (24) as 


4 LU M2 RR’ tt} ame 26 4 UM RSR’ 3+ cos aig B— log r—3) 
2 2 - y2 


U M?RR’ cos 26+4U(4M*—1)RR' logr+ 3ay?N+O(t*), (26) 


noticing that ¢d. = O(t?), as was assumed above. 
Finally, equation (8) can be written 


20304, 
; 3\ / 322 4UB = B\ /4U 
— { -*)| : 1B es = cos 26}cos0—r-¥{ r+ ‘| 7S hn 20)sin 8+ 
ye r° r r r2 
+O(t) 
Os sf s 36 
= 4p898" _ p77 p28". gy2p ~~ +0(1). (27) 
yr? i r3 
\ solution is 
08 G os @ . 208 ¢ s 36 s 
Mad, = $B?2" _ UB? 42 B= aI pce + QM" + 018). 
r? re a r r3 
(28) 


Here P and Q are arbitrary. Inserting (28) in the boundary condition 


8a) we deduce that 


ae 5 B 3U B =~ 3P Q 
— I os b— O08 pea OS 7+ 1] 2 y S< Diet y, Ss; —— 208 
q Ui 6 Fe é Ra °° 6+-4 pe 36 pa 008 30 R°° tf] 
- O(t?), (29) 
3M?—2 ; rs 
or g =" —auR, P= UR. (30) 
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Hence 
cos 8 —6 cos Y 


4 cos 30 
¢, = U M*R*—— —+2U M?R* or ao - 


(13 _—— ‘2)cos 6- e—3 3. M*cos 30 
r 


+i3U BR? 





—, (31) 
neglecting O(t). 
Hence, if O(t®) be neglected, we have on the body {r = R(x)}, by (3), 
Q; = 2U[4U M?R’ cos 6+4U M?R' (cos 36—6 cos 6) + 
+3U R’'{(13.M?—2)cos 6—3M?cos 36}|—2U?R’ cos 6+ 
+2U R'| —U M? cos 6—1U M*(cos 36—6 cos 6) — 
—,U{(13M?—2)cos 0—3M?cos 36}— 3U cos 0|— 
—2.2UR{AUM?RR' sin 26—2U M?RR'(—4)sin 20— 
—§U M*RR’sin 26} R-*sin 0 
= U*R’'{(4AM?—§8)cos 9A—Y.M? cos 36}. (32) 
We now apply equations (11), (20), (21), (22), and (32) to obtain that, 
at a point on the body, 
C, = {-—2U-(K'+A’' log R)— R”}+-¢ cos 0(—4R’)+-?(—1+-2 cos 20)+ 
+)|_R'{(—44.M2-+-§)cos 0+41.M?cos 36}+-4.M? 4R’(1—2 cos 26)cos 6], 
. (33) 
if O(s*), O(yFt8), O(?t?), O(yt?), and O(t*) are neglected. The resulting 
value of D,, by (12), is ; 
Dy+ po U22x | R(x) R'(x)(—Y) de, (34) 


0 
where D,/7p,) U? is 


Zz 


l 
— | BR (20- “S er : Ah dy +204 log E+) dx 








I x 
) 
=iR Ro ft l—y)A'(y) dy — p46 yar | Yow wh wha 
- i ie (Roe | de 
dx 2 
0 
iff 
= Ti [| A'(x)A'(y)log — dady — = R, Ro | A’ y)log — dy+ 


So 
So 
o 


D) 
+ R? R@log —., (35) 
3 
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where Ry, Ry are the values of R(l), R’(l). Expression (35) is of course 
the usual one for the drag at zero yaw, reducing to the first term only 
when either Ry or Ro is zero. By (34) we see that the addition to D, due 
to yaw is —}p, U*7 R§y* to a first approximation. 

The value of L,, by (13), is 


l l 
Log Uar\ ob R(x){4R'(x)} dx+y8 | — 
0 0 


= 4p) Ua R24 2-+(ZM2—$)y}. (36) 


Here terms O(/t*) are neglected. The true lift L is (36) multiplied by 
(1—4?) plus terms (from the drag) of orders which we are neglecting. 
The correction to the lift coefficient C, = 2% for thin bodies at higher 
angles of yaw therefore is ja. The percentage correction is 707?; the 
most important thing to notice about this is that it is small—e.g. for 
M = 2, & = 5°, it is 2-7 per cent. Thus we see that for thin bodies of 
revolution the lift coefficient is closely a linear function of angle of 
incidence. 

The true drag D, by (15), is D,+D,+4p,) U?7 R224? to the orders we 
are contemplating; thus we can write Cp = Cp,+¥? for thin bodies with 
a flat base: for bodies pointed at both ends, however, Cp is not changed 
to order ?. 

The moment is easily seen by (16) to be 

: 
Loy U2af4yb+ (34M2—8)f} | x R(x)R' (x) da, (37) 
0 
if terms in yt* and y¥*t* be neglected. The percentage ring in Cy, for 
7M? 


thin bodies at higher angles of yaw is therefore =~. This also is 
) 


small; for M = 2, % = 5°, it is 3-05 per cent. 


4. Second approximation to 4, 
Equation (6) can be written 


2a3(V*d, —r-2d,) (U 3) (2 U a —25)+<4(-" 4 ae 


? r r3 r ry? r4 








| uf A'_,AB' <=) + 0() 
r rT 





8A2B 12Uz 2+4U y. 8U2A’ 
84°B 12UA*+4UA'B _ +.0(0). (38) 


r°? ra 
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This can be solved by assuming that 


& = -+( +- a Frlogr) +-O( 5). (39) 


if C= A®B/2a2=}UM?R*R?, D = UM%(4R?R2+ RR’) 


and F = (2M?—1)A’. We obtain G from the boundary condition (6a) 
which takes the form 


_ BO G..ie8, 2. mr, B’ 
l Ri Rit D R +H+F(1+log R)=R Rt O(t*), (40) 
giving 
G- = 4. D(1—log R)-+ ER?-+ FRX 1+log R)-—-B’RR'. (41) 


The determination of H is more difficult. If equation (6) be linearized, 
it takes the form 


2a3(V%4,—1-24,) = 4U%Pgry-+ 2U azn, 
2 c a2 l 
or ie: es (42) 


ore or OF ac? 

This is the equation satisfied by the coefficient of ¢? if 4, is expanded in 
ascending powers of ??. The ¢, on the right may be replaced by the 
linearized ¢, which is — K,(apr)A(x). Now equation (42) must be satisfied 
by those terms in (39) whose coefficients are O(t?), i.e. Br-!+ Er-+ Fr logr, 
together with terms involving higher powers of r. Thus to determine £ 
we must find the necessary relation between the coefficients of r-! and r 
in the expansion for small 7 of a solution of (42) which satisfies condition 
‘CO’ at infinity. To do this we write (42) in Heaviside form (remembering 
that ¢, = 0¢,/éx = A(x) = 0 at x = 0): 


ard li . 
Or (se + at = 2M*?ap*K,(apr)A (x). (43) 
or? r Or re 

(Watson’s notation K,(z) —K({(z) is used.) A particular integral satis- 


fying condition ‘C’ (i.e. behaving like e-*”" not e%”" at infinity) is 
dp; — K,(apr) ( {2M?ap*K,(aps)A(x)}sI,(aps) ds— 
0 
I,(apr) [ (2M*ap?Ky(aps)A(x)}sKy(aps)ds, (44) 


since the Wronskian of J, and K, is —1/r. Equation (44) can be rewritten 


d; = Ky vpr) [ 8K,(aps)L,(aps) ds+L,(apr) ( sKY(aps)ds|| 2ap?M?2A (x) |. 
0 r (45) 
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The expression in curly brackets is, for small r, equal to 





j * 3 s))* - . 
— | ( Japs ds+ Japr| —[sK,(aps) =P") | sKi(aps)ds + 
xpr J \aps}” . | Pp Ie . 
+O(r?logr-1). (46) . 
But ( sk% vps) ds 1/2 x, so (46) is 
0 
9 
i pr ; ;| log ——y]— : | 4 O(r3log?r-), (47) 
dap “ \ 2p? ypr 2a%p?| 


and (45) is 
d; pr(log = 7) M*A(2)+-O(r* logs yt). (48) 


To (48) can be added (in order to make the coefficient of r-! correct) 
the expression 2«K,(apr)A(x), which is finite at infinity and for small r is 
asymptotic to 2A(x)/pr = B(x)/r. The result is 


br 2 } + tapr(log + y—4)}-+M*prlog "2 +7) |4(a, os 
xpr 2 Z 


in which the coefficient of 1/r is B and that of r is 


E (2M Plog ‘* + ”) a sp A(o) 


4a?A’(2x)- eMe—1)= { A W)log 





dy. (50) 
0 
We have now determined all the coefficients in (39) and can proceed 
to obtain the coefficient of %cos@ in C,. By (11), (3), and (6a), this, at 
a point on the body, is 


| Bo O' 4 D'logR. ,, 
— 2(T + + o°g ?) iter == © a y ‘“ y 
pa{2(U+K' +4’ log B+ + 5+— se" + ER+E Rlog R) + 
Le ee enema A\(6UA 2AB 
ry 4. 4M2U-4[2U(K'+A’ log R) 4 4\[2"* _ 4") 4 oes 
R Ri"? (2U °6 \+Fa)( Rk wy oO") 
aif .C GS pies, wes » , 
ae 1 - R 7 R3 -— R a R —— E R+ F R log R) = 


4x? R’ 


7 (K'+-A’' log R)+-2(M?—2)R3+0(8). (51) 
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The coefficient of 4% in L,/zp, U? is therefore (by (13)) 


l l l l l 
U-1 ({ 2 B dee-+ mete | G'de-+ | D'log Rde+ | BR de+ 
‘0 0 0 0 


1 


9 








l ‘ 
4 | F’ log Rar) — 
0 0 
l 
— R244M2R3 R24 { M2 RR’3 dx— 3M? R2 R2-+ M2(4R3 Re2-+ RB Ri)+ 


0 
l 


2 R2 df a 
| sar ()-+-(2i0*—1)— | A’ (ejlog—_—d 
+75 ( era’) +QM—1)5 | 4’) cea’) + 
0 
+ R§(1+2 log Ry)(2M?—1)(RP+ Ry Ro)—2R§ Re— 
1 


—M? | (4RR3+ RRR") dxe— 


A w(- bay? A’(x)+(2M?2—1 es [ A'(y)log —" i) da— 
‘ 4 sa . 


) 


U2 


Stim, Oe 


l 
~(2M2—1) | (R-+_ RR")R(1-+2log R)R’' dx— 


0 
| : / 1 - 
ona [ Aca [ log wie A'(y)dy+A'(x)log 4) dx+ 
1 x—y 2 


ax y 
0 


“ 
0 


1 
— M7?) ) | RR®° dx 
0 





Ll 
8VM2_9 a 
= R24 3M — J }ee A’'(x)A'(y)log — at. dady— 
U2 ~2la—y 
(5 2.9 J ’ 
pon 0.9) Fa Re l A ‘(y)log Ry dy+ 
U , 2(/—y) 


0 


2(2M2—1) 





aR, 
a 2 A’ 
+ il ma | (y)log Ay) dy+ 


+ (3M?—2)R2 Ri2+2M*R3 R’. (52) 


‘+A’ log R)dz+(2—M?) | RR? dx 
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The coefficient of % in the true lift (14) differs from that in LZ, by 
D,-+ Doo, Where Dy is given by (35) and D,y is the base drag at zero yaw: 


thus it isT 


Lt 
Tt? Dé 9 - R , ’ 
mol | Re 3a* | | log a d{ R(x) R'(x)} d{ Ry) R'(y)}— 
, 00 ; 
"oR 
602R,. RX log — ae 5 Rt ' ee 
a? Ry Ry | 85) URWIR (y)}4 


0 
l 


2(2M?2— 1) R al log . a{RW)R)}) + 





dl\ J ° 2(1—y) 
‘0 
| (3M?2—2)R2 R2+2M? RS B} ~De (53) 


This gives the correction to the formula L = zp, U?.R2 due to terms of 
order t*; earlier we found the correction due to terms of order ¥*t* to 
be ia*b zp, UR. 

For a body pointed at both ends we have simply 

oe 3a2Dp, (54) 
dis 
by letting R, > 0 in (53) and comparing with (35). 

If the pseudo-lift L, on the portion x <I’ of the body is L(l’)—given 
by dzp, U? times (52) with l’ written for 1, R(l’) for R, and R’(l’) for 
R,—then by (16) and (13) the moment M is 

1 

f LI 


l 
67 + Re) R'(«)}— da L,(l+R, R)— | L(x)(1+- R'2-+RR") dx 
, 0 
l 


l 


L, l+-udrp, U? RB Ry L(x) dx—ayrp, U* [ R?(R?+ RR") dx+ 
0 0 

+O(t®)-+O(2), (55) 

an expression of which, for brevity, further simplification is omitted. 

The term in #2 was given in §3. As an example we consider a cone of 

unit height and base diameter ¢t. Then 1 = 1, R(x) = }tx. We find that 
i 

jp, Unt? 


M2 ¥ 3M2~— ” 
a (2 ae P—Co, i+ FaX2-+- OC!) + OW), 
(56) 


+ Here and in (52) d/dl does not operate on Ry. 
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and 


“a M 
M™ 1p, U2 i nt? 1 
2Po qm 


2 
(1 + LM*log + pre) i 


—8 
9 fF + O(pt*)+ O(t?). (57) 
Tables 1 and 2 give the values of dC,/ds at 4 = 0 (omitting the base drag 
term), and of dC,,/d/ at 4 = 0, respectively, for M = 1-5, 2, 2-5, and 3 
and for cones of semi-angle « = 
may, of course, be expected to 


5°, 10°, and 15°. The values for smaller « 
be more accurate. 


TABLE ] Sm TABLE 2 








we 
(Cf. 2-000 at « = 0°.) (Cf. 1-333 at « = 0°.) 
M e= 5 €= 10 €= 15 e=5° | e=10 €e= 15 
“s 1°894 1-711 1°526 1298 | 1°257 1°247 
20 1°887 1758 | I°741 1289 | 1270 | 1°348 
7 758 | 74 | | 1°34 
2°5 1°881 | 1°831 | 2°045 1°282) | 1°308 1°524 

5 : : | ‘ : 

370 1'879 1°945 2°473 1°279 | 1375 | 1°790 








It may be noted that our use of the exact equation does affect the terms 
in ¥t? in CL, since if the linearized equation be used (1), then dC, /dy isa 
function of at only. The values of dC; /dyb given by the present theory and 
by Tsien (1) are shown in Tables 3 and 4 respectively. (The good agreement 
at M = 1-5 is seen to be merely due to the intersection near this point 


of widely diverging curves.) 











la a - @ r 
TABLE 3 [TABLE 4 
Msie=s5 ¢= 10 €= 15 e=65 ¢€= 10 «= 15° 
rs | 1°932 1°824 1°727 1°944 | 1832 | 1°700 
2°0 I°9QI9 1°843 1°879 1°887 | 1'692—s| 1°492 
2°5 1°908 1°899 2°143 1828} 1°565 1°324 
| au 595 : 

30 1'903 | 2°000 2°541 1°768 I*451 1°186 








5. Conclusion 
The pressure coefficient on the body (by (33)) is 
O(t2) + O(t) + O(p2). 

If, as seems likely, this is true even when the presence of the bow shock 
is taken into account, we may deduce that the strength of the shock is 
at most of this order.t Since the variation in entropy behind a shock 
is known to be of the order of the cube of its strength, we may conclude 


+ In a later paper it will be proved that when 
actually O(t*). 


0 the strength of the bow shock is 
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that the error in using the adiabatic equations of motion behind the shock 
is at most 
O(t®) +- O(et®) +- O(xp*t*) +- O(F #3) + O(4t?) + O(>t) + OY) 


in the pressure coefficient. We have neglected terms of these orders and 


oi) 
, may therefore deduce that our expressions of order ft? and 4° in C,, and 
drag | Cy additional to the known term of order y are correct. The results were 
ad 3 given, respectively, at the end of §4 and at the end of §3. 
ler € 
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SUPERSONIC FLOW PAST SLENDER BODIES OF 
REVOLUTION THE SLOPE OF WHOSE MERIDIAN 
SECTION IS DISCONTINUOUS 


By M. J. LIGHTHILL (Department of Mathematics, 
The University, Manchester) 
[Received 10 November 1947] 


SUMMARY 

The theory of supersonic flow around slender bodies of revolution, yawed or 
unyawed, with pointed or open bows, based on the linearized equation, is extended 
to the case when the meridian section of the outer surface has discontinuities in 
slope. Expressions for the pressure distribution on the surface are obtained. It is 
found that the drag coefficient is no longer independent of Mach number, and tends 
to zero more slowly than the square of the thickness of the body. The large pressure 
change behind a discontinuity is made up remarkably rapidly. The first approxi- 
mation to the lift coefficient is unchanged. 


1. Introduction 
THE study of supersonic flow past bodies of revolution was initiated in 
1931 by von Karman (1), who used the ‘linearized equation’ of com- 
pressible flow. He advocated the exact solution of this equation under 
the exact boundary conditions—this involves the numerical solution of an 
integral equation—but also gave an approximate solution for the case of 
a body pointed at both ends. In 1945 the author (2) removed certain 
anomalies in the use of the linearized equation in this problem, and 
established von Karman’s approximate solution as a first approximation 
both for bodies pointed at both ends and for shells (in the region ahead 
of the flat base). Shapes with minimum drag were given and a theory 
was also evolved giving the flow outside any duct of shape approximating 
to a circular cylinder. The theory given for the flow inside the duct was 
incorrect. Broderick (3) (1947) has obtained a second approximation to 
the flow past pointed bodies and shells, using the exact equation of 
adiabatic flow: terms due to use of the more exact equation are shown 
to be important and it is inferred that, where the linearized equation is 
employed, not much is gained from proceeding beyond a first approximate 
solution. Meanwhile in 1938 Tsien (4) applied von Kaérman’s method to 
the yawed flow past a body of revolution. In a recent paper the author (5) 
obtained correct second approximations in this problem. 

In all these approximate solutions for thin bodies it is necessary to 
assume that the slope of the meridian section be continuous and that its 
curvature be bounded. When the slope is discontinuous, the steps by 
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| which the solution is reached are invalid, the solution is wrong, and, in 
lfact, the expression which it gives for the drag of the body is infinite. 
Nor can these anomalies be removed by smoothing out the discon- 
tinuities by small curved portions whose curvature is allowed to increase 
ndefinitely. 

In this paper the theory of supersonic flow (both symmetrical and 
awed) past thin bodies of revolution is modified to allow for discon- 
tinuities in slope. The theory is presented in such a general form that 
talso applies to flow outside slender ducts whose thickness is small com- 
ared with their length (but whose radius is not now necessarily approxi- 








jeqd | mately constant). Approximations to the pressure distribution, lift, drag, 
;in |and moment are found for the flow past any slender body of revolution 
‘is | the outer boundary of whose meridian section is made up of a finite 
1s ; 
: umber of analytic arcs. 
ure . 
«. | 
)XI- 
2. Lemmas obtained by Heaviside calculus 
By ref. (2), §1, or ref. (5), §3, an approximate disturbance-velocity 
potential in the symmetrical flow will have the form 
in 
r—-ar 
m ‘ ft (y)dy , 
dy = a = Ko(xpr) f(z), (1) 
el : \ (x y)* — i 
an ” . 
ol vhere a = ,/(/?—1), cy is any positive constant, p-! is | dx (in Heaviside 
4in a ; 0 
all notation) and Ky is the Bessel function as defined by Watson. For small r, 
on | |!) 18 approximately equal to 
ad * Din 
a 2(x—1 
me flog(2/apr)—y} f(x) = log {2 WN afiy) (2) 
ry | or J 
ng ° 
as y being Euler’s constant), by the Heaviside representation of a Faltung 
to ntegral and the fact that loga = —logp—y. (We assume f(x) = 0 for 
of | «<¢y.) If f(x) is continuous, and zero at x = Co, the derivatives of (1) are 
wn cee 
is Od, y f(y) dy . . 
| Sf a ay — PKlapr)f(e) (3) 
ite O2 J Ue —Yy)P—aer*s 
to Co 
5) ads r—ar 
; Og l ff («a—y)f'(y)da : x 
| nd 9 _— | DEY — —apK(apr)fiz)~ 22) (a) 
pn ee TF) ey aPrh / 
to Cc 
Its asr—+>0. (On Watson’s definition K,(z) = —K,(z) ~ 27! as z > 0.) 


by In ref. (2), §5, a function U(x) is defined (and tabulated) such that 
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(equation (78)), if in (3) and (4) r is given the constant value R, then the 
relation 





or 


we deduce that the Heaviside representation of U(x) must be Ko(p)/K,(p). 
(This fact was first pointed out in an unpublished Admiralty Computing 
Service report. A complete theory of the flow outside and inside tubes 
of approximately cylindrical shape, based on this approach, is given by 
Ward, ref. (6).) Hence, for large x, | 
x . ; 

U,(z) = | U(y)dy = Kol) - log(2 PY—Y O( p* log p-1) 

pK,(p) p(p~) 


0 
= log(2a)+ O(a-* log x) (7) 
and U(x) ~ z (8) 
x 


For the yawed flow at angle 4 to the axis, the disturbance-velocity 
potential can be written $9(7, x) +%4,(r, x)cos 6+ O(x?) and an approximate 
equation for ¢, is 





a2 a2 
ed, | oi si vse ia = 


or rr Ox? i ) 


or? 
a solution of which is 
r—-ar j 
l (x—y)g(y) dy 
=> | 3 


= = aK,(apr)g(x). 10 
1 r J{(x—y)?- x22} 0% (apr )g(2) ( ) 


ry 
Co 


The derivatives of ¢,, if g(x) is continuous and zero at x = Co, are 


x xT 








a a | (e—w)g y) dy ee | ofa) (a) | 
Ox re J A(xz—y)*—a*r*} r 
(as r > 0) and 
“P (w—y)°g’ - 
as | (7—y)"'9 (y) dy = ap Ky (aprig(x) ~ —— | g(y)dy. 
or re J i(z—y)*—a*r) aa 
Co Co (12) 
We shall need to solve the integral equation 
z—aR 
f = (w~—y)?g'(y) dy eK’ 
s—¢—aRk — . vi! _ — —g2nR?K' (apR)g(a (13) 
_ V{(e—y)?—0? R?} sical 





(for x > 
) solutior 


Ob it | u (27=Y\ q(2bol"-) (5) | 
or té‘i‘i ‘ xR or may 
0 


subsists. Since, by (3) and (4), the relation must be 
py oe Ko( xp R) Oho (6) 
ox VK (apR) er’ 


where 


as p > 


The fu 
ref. 6): 


| 3. The 
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for « > c+aR), and find é¢,/éx, when r = R, for the g(x) which is the 


solution. Clearly we have 


od, : : p-te-c—ak 1 ,/x—c—aR 
— I »R) > - = F — I, 14 
oar PANOP") 2 RK (pk) PR xR i 
where F (x) K,(p)/pKy(p) = 14+ O(p*log p-) (15) 
is p > YU. Thus 
F(x) >lasxa—-oo, and [ {1— F(x)'dx = 0. (16) 


0 


The function F(a) was first used by Ward, who denotes it by L(x) (see 
ref. 6): a table is given in § 4. 


3, The symmetrical flow 

Using cylindrical polar coordinates (r, 4, x), let r = R(a) be the equation 
f the outer boundary of the meridian section of the body of revolution 
onsidered; let it extend from x =a, to x =a, and have maximum 


thickness t. Let R(x) be continuous in a, < x < a, and analytic in each 


interval ag <X% <@,,0, <2 <Qy,...,@,_. <% <a,. Putt 
R'(a,+0)— R'(a,;—0) = 5; 

interpreting this for 7 0 and 7 n by the convention that R’(x) = 0 

for <a, and x>a,); put also R(a;)= R;, a,—aR;=c; (where 
/(M*?—1), M being the Mach number of the undisturbed stream), 

S(x) = 7 R(x) (this is the cross-sectional area at distance x from the 
rigin); then we must have S’(a,;+0)—S’(a,—0) = 2b; R,;. 

When there is no discontinuity the velocity potential is given (to a first 
approximation) by equation (1) with f(x) S'(a)/27. We take as our 
corrected value of f(x) the expression 


r 


ae <¢ o 
f(x) = —z S"(y)dy— > f(x). (17) 
“a7 . 1=0 
vhere f,(z) is zero for 2 < c; and continuous for x >c,;. The integral here 


which can be written p-!S"(x)) is not a Stieltjes integral with S’(y)dy 
written for dS’(y)—we shall always write Stieltjes integrals as such 
throughout this paper—but is an ordinary Lebesgue integral which 


i 


ignores the points y = a; where the integrand is not defined. Thus in the 


} 
interval a, v a,,, it is S’(x) S 2b; R;. 
. i=0 
[he boundary condition to our order of approximation is that 
C>,/er = R'(x) when r R(x). Using the approximation (4) for the part 


Dashes denote differentiation with respect to x. 
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of é¢,/ér due to the first term in f(x) we can write the condition (in the 





interval a; <4 <4,,,) as 
1 j 
No atte. ‘(x)— ¥ 2nd, R, 
R'(x) = oh R(x)R' (zx) p: nb, 4 


(a—y) fily) dy 
BD tz Key Rey 


Equation (18) can be satisfied for all 7 if, for all i, f;(7) satisfies 
xr—aR(x) 


(x—y)fi(y) dy 
b, Rk, = Ce a SO J 9) 
—— | \ {(x—y)®—a? R*(x)} i" 


Ci 





The term P, due to f;(x) in é¢)/éx is then given by (3) with —/f; written 


for f and R(x) written for r. When zx is near a;, so that R(x) = R,, an} 
approximate solution, by (5), is 


P= — 210 Pe 20 
(making the convention that U(x) = 0 for x < 0). Rete x is not near 
a; equation (19) becomes, by (4), approximately f,(7) = 6; R;; and in 
r—aR(x) 
; pri (21) 


i a2 R? 
NV {(c—y)?— R (a °)} 

Cs 

the greater part of the integral comes from a small region near y = ¢, in 


which f,(y) rises rapidly from 0 to b; R;; so that we can write 


U 
OR; - 
—  2—a, (4) 
Since (20) also satisfies this approximate equality (by (8)) when z is not 
near @;, (20) can be taken as a good approximation to P, for all 2. 
We deduce that the pressure distribution is given by 





D a ¢ 
vs or 


kp, V? 6a 
x—-aR(x) y’ 
on [ Shy) dy +> 26; 1 (241) _ pryx) (23) 
7 J Ni {(a—y)? — oa? R?(: yx 
ao = 
on the body, approximating as in ref. (2). 
It should be noticed that (to our degree of approximation) as soon as 
each summed term in (23) attains its asymptotic form 2b; R;/(~—a,), it 
becomes the term due to the discontinuity at y =a, in the Stieltjes 


= PBs atts) 





integral 


xr—aR(x) 


1 dS" ( (y) 
7m J {(e—y)?—o? R%X(x)}? 


0 


(24) 
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which is the natural extension of the form which applies when discon- 
tinuities are absent. Thus we may say that each discontinuity (including 
the front if S’(x) #0 there) affects appreciably only the region imme- 
diately behind it: there is a sudden change of pressure (a drop, in the 
case of an ‘expansive’ discontinuity); but in a distance of the order of the 
thickness ¢ of the body the pressure is made up to what it would have 
been had the shape been continuous. That this effect can be appreciable 
is indicated by the fact that the sudden change in C,, of 2b,/a is of order t, 
while the other terms in (23) are O(t#logt-1) at most. 
Writing D for the drag divided by 3p,V?, we have 
D = | C,S'(x)dz, (25) 
0 
where C,, is given by (23). With U,(x) defined as in (7) the term 


x 


- _,[x—a,; 
U- ‘) S’(x) dx 26 
| ( ze) (x) (26) 
U 
i r—a; 
becomes ry | U,i—— 968 (2). (27 
- ; ( xR, ) (x) ) 
Now by equation (7) we can write |U,(x)—log(2x)| < Aja~loga for 
a > A,, while for x < A, we have |U,(x)| < A,, say. Hence (27) can be 


written 


- x—a,\ - P x3 R3 xr—a, 
—_ ool? : S (7)+O ade. ST eR: ‘ V(a)| + 
R, | log (27 )¢ (x) | los( ; Jia (x)|+4 


i . (x—a;,)? x Lt; 
a; aj+ A ,aR; 
a+ A, aR; 
x—aa,| : 
ais vR, (1+ log2——“! | \d8’(x)|], (28) 
my rR; 
and the O-term is O(t*log t-). 
Again we have, by (2), 
YR(xr) Z 
i Ss” | y) dy [ 2(x — y) - 
pons — | los—_~" dS" (y), 29 
4 \{(x—y)®— a? R2(x)} > aR(a) (y) (29) 
0 


correct to order ¢? except in small regions to the right of each discontinuity 
of length O(t), where the discrepancy will not affect the integral (25) to 
order ##. Multiplying (29) by S’(x) and integrating from 0 to 00, we obtain 


a 
9 


0 zx 
S"(ax)S8" (x)log —— da + | S'(e) dx { loge —y) a8"ty), (30) 
x R(a) . 
0 


. 


0 0 











96 M. J. LIGHTHILL 


and the second term can be rewritten 








( dS" (y) ; S'(x)log(x—y) dx = —| S”(y) a( [ 8'(e)log(—y) de) 
0 y 0 y 
— | S”(y) dy | log(x—y) dS’ (x). (31) 
0 y 


Collecting terms we deduce that 


O @ oO i.2) 
D . [ $'@)s"@)log- - -dx—* [ s*y)ay [ tos(e—y) d8"(x)— 
7 ~ wR(a) T | 7 
0 0 y 
2b; R; | log "hs dS'(x)— | S'(x) R(x) dx. (32) 
i=o rR; 
aj 0 


The first and fourth terms combine (since S’(x%) = 27 R(x) R'(x)) to give 


n 


n 7 
‘J d | ’2(a)log +5 _| dx = ss > (los a) S%as—0) —S'*(a;+0)}, 


on | dx\ (a)J an £4 i 


(33) 
and this combines with the part of the third term 


nN. oo 9 1 n 2 " - “ 
_ 2 2b; R, | loge 48'(e) aa > (los alts (a,+0)—S'(a;—0)}S'(a;--0) 
i=0 de cl i=0 : 


(34) 
to give 
l n 9 9 
= SQ’ os —s, — = 92 oO ae . 5 
a > (oe as (a,+0)—S'(a,;—0)} 2a S b? R? log Ry (35) 
i=0 i=0 
The remainder of (32) can be written 
am | S"(y)dy | S"(x)log(a—y)da— > 2b; R; [ s"wilozia, —y)dy— 
T i 7 
0 y : 0 
n 
— 2, 2B ‘\, Js r)log(a—a,;)da+ 2, 27b; R; log(a;—a;) 
Gn an n An 
= : | [ S"(x)S8"(y)log L _ dedy 4. S 2b; R; [ S”(x)log - dx + 
ye L—Y — J L—A;| 
ag ap a ao 
— ] 
YS 47b,b, R; R,log- 36) 
he de ial "*ia—a,l” ( 


and D is given by (35) plus (36). 
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The form of (36) is not surprising, being a natural extension of the form 
of D when discontinuities are absent. For consider the formal Stieltjes 





integral 
aes i l 
31) - | a8’ (2) | log ———dS'(y) 
2a J t—Y| 
0 0 
“ x n @ 
} ] F - - l =~ J i ve 
——_ | dS'(x) | S”(y)log —dy + > — | 27b; R; log ——— dS’ (zx) 
2a J x—y| ~ a+ Dz x—a,;| 
0 0 = 0 
1 -" y 1 m - 1 
= — | S” (a) dx | S” (y)log ——— dy +2 > b, R; | S’(x)log ——— dx + 
2a x—y " |x—a,; 
29 os os ; “— ao 
Ju) n n . 4 
+ > b, R; > 2nb; R;log———.__ (37) 
— 5 
, i=0 j=0 si e 


If the infinite terms with i = j in the latter double sum are all discarded, 
(37) becomes identical with (36), which can therefore be appropriately 
| termed the ‘finite part’ of the double Stieltjes integral. Using an asterisk 
33) | to denote finite part we may then write 


x x 





nN 
l*r an Sie 2 
D=— | | log ——d 8’ (x) dS’ (y) +2 D3 b? R?log——-, (38) 
a x—y é ak; 
UY) | 0 0 iis 
where the first term is defined by expression (36) and has properties similar 
(34) to those of D when discontinuities are absent. The second term, however, 
is quite new and makes two important changes: 
,.. ) (i) D is now dependent on Mach number, decreasing as the latter 
(Oo) . 
m InCTeases ; 
(ii) D now contains terms in f*logt-! as well as ¢4; the former are of 
course greater for small tf. 
As an example consider the flow past a double cone of thickness ¢, 
lefined by the equations 
ta 0O<2z<} ; 
R(x) =! ( 2)\ (39) 
| \t(l—x) (<x <1) 
There is one discontinuity, a, 4, where b, = —2t, R, = }#t. We have 
S"(x) = 2nt? for 0 <a <1. Equation (23) becomes 
1 (2 cosh-1(1/at)—1 (0 <2 < }) 
' (7 — 
p 
| iol x 4,,(2e—1 (40) 
t | 2 cosh 1 J-1-<1 (4 < zt < 1) 
™ xt(1—ar), xl ott > 
(50) , 


and equation (38) gives 


Up = 
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C)/t?, and values of C,,/#? as a function of x, are given for five values of 
at in Table I. The relative drop in pressure at x = } is seen to be very 























TABLE I 
| 

xt O° ee 0-5 | 050 | 0°55 | 060 | 0°65 | 0-70 | 0°75 | o-80 | 0°85 | 090 | 0°95 Cpt 2 

— = sane | a in 
bi 62 — |—73°6| —29°0| —13°5| —6-4 | —2"4 | +03 | +23 | +40 | +5°6 | +7°7 || 32°84 

pe 
oS fe) + 5*0 ~35°0| —20°2| —12°1) —7-0 | —3°7 | —1-2 | +0°7 | +2°5 | +42 | +62 || 14-42 
OrI5 +4°2 | —22°5] —15°0] —10°0| —6°5 | —3°8 | —1°7 | +0°0 | +17 | +3°3 | +5°4 11°18 
0°20 +36 — 16-4} —11-8) — 8-4] —5°8 | —3°6 | —1°8 | —o-3 | +1°2 | +2°8 48 | 8-88 

> | i] | 

0°25 +3°1 Ea = ig) — Fal oF. | Ss | HPS | HO eo | 4-0°5 #5 || 7°09 











great especially for small at: the recovery is also very rapid. The maximum 
of dC,,/dx (which occurs just behind the discontinuity, after which it falls 
rapidly) is approximately 4/«?, since U’(0) = —4. When WM is near unity 
this may cause boundary layer separation or at least transition to turbu- 
lence: experiments with double-wedge aerofoils have not indicated separa- 
tion as likely when the angle turned through expansively is small. In 
ref. (2), §3, it is shown that for a parabolic shape the quantity C/p/t? is 10-67, 
for two other smooth shapes 11-64 and 12-5, and that the minimum possible 
value for a smooth body symmetrical about a central plane is 9-87. Thus 
we see that for a very thin body the discontinuity increases the drag, but 
that for thicker bodies this is not so. In fact we would be tempted to 
infer that for af = 0-20 and 0-25 the discontinuity lessens the drag, were 
it not probable that terms of higher order than ¢? in Cp become important 
for these values. 

We observe, in passing, the import of our results in the theory of 
cylindrical sound-waves produced by the expansion and contraction of 
a small cylinder: the mathematics is of course identical in the two theories. 
When the circumferential velocity is changed discontinuously by an 
amount V, the surface pressure changes by the relatively large amount 
pVa; but this change is largely made up in an interval of time O(R/a), 
where R is the radius of the cylinder: this interval is very short. Again, 
our expression for the drag of a given body of revolution gives the work 
done in expanding and contracting the small cylinder in a given manner. 


4. The yawed flow 

We first treat the yawed flow on the basis of equation (9), considering 
afterwards how far the use of a more exact equation will alter our results. 
When discontinuities are absent we have equation (10) with g(a) = S’(x)/7- 
We take as our corrected value of g(x) the expression 


gx) =~ | S*y)dy+ ¥ gia), (42) 
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where each g;(x) is zero for x < c; and continuous for x >c;. Our boundary 
condition to the first order (5) is that @¢,/ér = —1 whenr = R(x). Using 
the approximation (11) for the part of @¢,/ér due to the first term in g(x) 
we can write the condition (in the interval a; <x <a,,,) as 














}_ p2 =— i = 5 ) 
| SEN" (x) a Lo ni R, (a —a;)} — 
; «2 aR(x) 2g “ \d 
S a—Yy)gily) dy ' 
J, (43) 
0 J ee > — a? R*(ax)) a)}? 


which can be satisfied for all if, for all 7, when x > a,, 
@— XR xr) 
[2] = | Wi (y)dy _ 
< V{(x—y)?- x2 R2(x r)}? 


“0 
Ci 


2b;R,(x—a,)- (44) 


the square brackets indicating that the term within them occurs only 
fori = 0. Hence by (14), when 2 is near a;, so that we may replace R(x) 
by R; approximately, the term q; in r@¢,/éx (on the body) which contains 


q, 1S 

: x R «r—a 
q; = 2b, R, F(t) 4. |= (= 45 
a“ se (“e)+ E “aR, (45) 


When is not near a;, however, equation (44) becomes, by (12), approxi- 


mately - 
26; R,(~—a,)+| Re] = ( g;(y) dy, g(x) = 2b; R;, (46) 
and we have , 
2—aR(2) 
q (Ya Ey _ x g(x) = 2, By. (47) 


J {(v—y)?—a? R*(2x)} 
Ci 

Since (45) also satisfies the approximate equality (47) when x is not near 

a;, by (16), we can take (45) as a good approximation to q; for alla. We 

then have (see ref. 5), to our degree of approximation, on the body, 


dC od 
—P — —2-1cos 0 
dis Cx 


z 
cos@j1 f¢ |, (x R x—a 
wis = Ss" dy +- Ss; 2 — 0}. 
malt | W)dy + D2 RFE =e) °F =. ) (48) 


v 








1=0 


ao 
The functions F(x), F’(a) are given in Table II, condensed from a fuller 
table in ref. (6). Since F(0) = 0, d¢ ',/&s is continuous everywhere: as each 
F attains its asymptotic value 1 (which it does rapidly, though over- 
shooting it at first) the first two terms in curly brackets in (48) become 
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simply S’(x)/7. The third term, which is only of importance near the 
front, where it predominates over the other two, is that occurring (by 
itself) in the theory of yawed flow past quasi-cylinders, due to Ward. 
(It is of course identically zero when Ry = 0.) 

Writing L for the lift divided by $p)V* we have (ref. 5), neglecting 
terms O(t*log t-), 

















dL r rd 
——_= — | cos0dé | R(x)—"dzx 
dis dis 
C7) 
an zx 
] y’ xr—a 
— 27 [ — | S’(y) eo mrs +—F *(———"}} da 
J |a a aR, 
ao Ao t=0 
an Ry 
y _— Qa _, (x—a 
—2 | S'(x)dx+4n > re b,R = do Bef PF’ a0) de, 
mi ak, 
ao aj ao (4s 
9) 
TaBL E II " — pnts 
ee: 7 By use of (15)—the precise form of (16)—the 
x Fe) | F@) integral from a; to a, is seen, by putting 
viene + 
lohie) 0*000 | I*00o0 j 
& == Gy -Y, 
o*4 0°359 0°796 itakpy 

8 638 “600 fa hl . : * 

- py pine to be O(t4log t-!). The last integral in (49) is 
I*2 0°543 0°425 ‘ 

1°6 0985 0°287 mee 

2°0 1'076 0176 “aRo 

2°4 1*I29 0°093 - 7 {4,—Ag 

2:8 1154 0°034 ah, [ F'(y)dy = ay F (" RR 

3°2 I'l59 —o0°005 3 teat | 

3°60 I°I52 | —0°029 

4°0 1138 -0°042 — | 5() 
a Bed Bee = aR,+O(#logt). (50) 
45 I*IOI —0'046 3 ° 4 = C 

52 1084 | —o1042 Hence neglecting terms O(ttlogt-1), (49) be- 
5°6 1:068 —0°037 comes 

6:0 1054 —0'03I 

6: aia | —o-o26€ : + ¢ _ ee 

ia) gua | soem 2{S(a,,)—S(ap)} + 27-Ry = 2S(a,). 

§ 1'034. | —0'020 g 
7°2 1027 | —o-016 (51) 
76 I‘o2I | —o'0I2 
8-0 I'017 | —or009 This is the ordinary result for a thin body, 








which is seen to be true independently of dis- 
continuities. (The other novelty is that we have established it for an open- 
ended body not necessarily of ‘quasi-cylindrical’ shape.) Though the 
discontinuities produce changes of pressure-distribution of the order we 
are considering, these integrate out to zero in the expression for the lift 
(as a result of equation (16)). 

Similarly the moment is unchanged. In fact, writing M for the moment 
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about an axis in the plane of the nose (perpendicular to the stream) in the 
sense tending to decrease the angle of yaw, divided by }p)V?, we have 
(see ref. 5), neglecting terms O(t*log t-*), 


on An 











: : 0 
dM -- cos 6dé@ | (x a.) Rix)" Pdx 
dis t dis 
0 ao 
a x | n ie | R 2 = 
= 27 (x a,)| S”(y) dy + > 2b, R; F ie -j- — F we dx 
e | a * i 0 ak; Oo aR, 
An An 
=2 | (v—a,)S'(x)da+4a > b, RB; [ (eas) t— Or) _y\ der 
P i=0 p | rR; | 
s) p — 
= [ (e—ay" (eae. 
a J ak, 
, (52) 
The integral from a; to a, in (52) can be written 
= 
| (a;—ayt+aR;y){F(y)—YaR,dy = O(Plog*t-*); (53) 
0 
and the last integral in (52) is 
adn— do adn— do 
rEg an—Ao 7 
YRS | yF'(y)dy RS Ly FY)— Ye — | {tFy)—'Ydy 
0 0 
O(log t A): (54) 
hence finally (neglecting terms O(tlog*t-*)) 
an an 
aM 4,f ” of tm —_ -_ 
dds = - (x A)& (x) dx = @ {S(a,,)—S(x)} da. (5: ) 
ag do 


Again the interest is that we have established (55) generally for thin 
bodies even with an open front. The aerodynamic centre is at a distance 


An 


( ( = alae (56) 
R S(a,,) 


ao 
behind the front (by (51) and (55)). 

We must point out here that, as was shown in ref. (5) (equation (42)), 
our equation (9) is not the true linearized equation for ¢,, which should 
contain a term 2M2e%¢,/éxér on the right-hand side. However, using the 
value of ¢, obtained in §3, and solving the equation as in ref. (5), §4, it 
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may be shown that additional terms in dC,/dy are all of order flog ¢-1 
which we are already neglecting. Similarly no change, except to this order, 
results from using a still more exact equation for ¢, as in equation (38) 
of (5). 

As an example, we consider the yawed flow past a double cone frustum 
of thickness t, defined by the equations 


gt(x+3) (O<a<}h 
R(x) = (HETI) O <2 <4) (57) 
\i(g—2) (b <2 <1) 
Here a, = 0, with b, = 3t, Ry = }, and a, = 3, with b, -t, R, = th. 
Equation (48) becomes 
dC, oe +-ot F (4a/at) + F’(42/at)}/(a+4) 
—a—Psec? = {¢ Pa ae 4 ~ 
dibs [ Qatar oat (4a:/at)-+ F’ (4x xt) — 4at F{(2a2—1) ot}]/(3—a)} 
(58) 


and this quantity is tabulated below for two values of af. The lift for 
this body satisfies dL/dy = at?/8, and the aerodynamic centre is at a 




















| | 
al x 0] o'05 | ovr | o2 | 0°32 o'4 075 | o755 | o6 o°7 o8 | og | r'o 
ess | 
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o'2 2°00 | 1°24 | 0°72 | 0°38 | 0°38 *39 | "40 0°07 | —or18 | —o-48 | —0°58 | —0°58 | —0'55 

















distance $ of the total length ahead of the front. (This is because, as the 
table shows, the pressures producing lift are concentrated at the front 
while those at the rear produce negative lift.) 
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THE FORMATION AND ENLARGEMENT OF A CIRCULAR 


HOLE IN A THIN PLASTIC SHEET 
By G. I. TAYLOR (Trinity College, Cambridge) 
[Received 15 December 1947] 


SUMMARY 

When a circular hole is made in a flat sheet by a conical-headed bullet or by 
outward radial pressure on its edge, the metal near the hole piles up into a thickened 
crater. The mechanics of this deformation is discussed. The interest of the problem 
lies in the fact that the complete strain-history of each element of the sheet has to 
be calculated. This is because the ratios of the principal stresses at each element 
of the sheet vary as the deformation proceeds, so that there is no relationship between 
the stress and total deformation but only between stress and strain increments 
occurring during a small expansion of the hole. 

If b is the radius of the hole at any time, the strain is found to be elastic at points 
where 7, the radius, is 3-646. In the annulus 2-216 < r < 3-646 the strain is plastic, 
but comparable in magnitude with the small elastic strain when r > 3-646. In the 
annulus 6 2-216 there is finite strain. At the edge of the hole the sheet has 
thickened to 2-61 times the thickness of the sheet. 

Experiments made with lead show that the symmetrical deformation contem- 
plated in this analysis does not occur ; but an alternative unsymmetrical deformation 
is produced which calculation shows to require less work, in the ratio 2-6 to 1-0, 


than the symmetrical mode. 


Introduction 
It is a comparatively easy matter to formulate simple laws of plasticity, 
but it is difficult to solve even the simplest special problems. Perhaps the 
most difficult feature in solving such problems is to connect the stress with 
the total deformation when the ratios of the principal stresses at each 
element of the material are not constant during the straining process. 
Most of the problems to which correct solutions have been given are 
concerned with plane strain of incompressible material, and in that case 
this difficulty does not arise. The chief interest in the solution of the 
special problem here discussed is that the ratios of the principal stresses 
change during the straining so that the whole history of the straining of 
each element has to be followed from the beginning to the end of the 
straining period 

In only one previously published case has this difficulty been faced. 
The present solution was obtained before that work was started, and as 
its authors pointed out, their recognition of the fundamental difficulty in 

+ R. Hill, E. H. Lee, S. J. Tupper, ‘The Theory of combined Plastic and Elastic 


deformation with particular reference to a thick tube under internal pressure’, Proc. Roy. 
Soc. A, 191 (1947 
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this kind of work was prompted by some unpublished work of mine which 
they had seen. The solution which follows is part of that work. The 
remainder will not be published as it has been superseded by the much 
more complete work of Hill, Lee, and Tupper. 


Enlargement of a hole in a sheet 

When a hole is enlarged in a sheet of metal by radial pressure, as it is 
when pierced by a pointed conical broach or a pointed bullet, the metal 
near the edge of the hole is thickened over an area extending out to radius, 
say, ‘2. Beyond this radius the metal suffers only small thickening, and 
outside some radius 7, it suffers only elastic strain. 

A simplified analysis illustrating the main features of this state of affairs 
was given in a private communication to the author by H. A. Bethe. With 
his permission a shortened version of this analysis is given below. 


Bethe’s analysis 
rT ¥ ° . 
The following assumptions were made: 


1. The sheet is thin and the radius of the hole is large compared with 
the thickness. The sheet is of infinite extent. 


bo 


Outside a radius r,, the sheet is in a state of elastic stress due to 
radial pressure at 7. 

3. The stress o, perpendicular to the sheet is zero. 

4. In a region 7, <r <1, the plastic strain is small so that the stress 
distribution adjusts itself till a flow condition 


og—o, = Y (1) 


is fulfilled. Here Y is the yield stress. 
5. Mohr’s condition of plastic flow is satisfied so that the maximum 
difference in absolute value between the principal stresses is equal to 
Y. It is further assumed that Y is constant for large as well as small 
strains. 
Using the equations of equilibrium and elasticity with assumptions (2) 
and (3), the following expressions for the stresses in the elastic region r > 1, 
are derived: 


radial stress o, = —A/r* 
tangential stress og = +A/r? } (r >74). (2) 
stress normal to sheet o, = 0 


Here A is a constant for any given state of stress but varies as the straining 
proceeds. 
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Using the equation of equilibrium and the plasticity condition (4) 


5) (r, >r >1.). (3) 


A _ 
— = }tY. 4 
e - ( ) 


Bethe pointed out that, if Mohr’s condition (5) is assumed, equations (1) 
and (3) are only applicable so long as gg is positive (oc, is necessarily 
negative), for it is only then that the maximum difference between the 
principal stresses is that between o, and og. When In(r/r,;) < —} the 
maximum difference is between o, and o,. The radius r, at which og = 0 
is therefore rs ret = 0-606r,, 


or ry 16575. 


Within the radius 7, the condition of maximum stress difference could 
still be satisfied if o7 = 0 and o, = —Y. The equilibrium equation 
applicable to sheets of uniform thickness, namely, 

do, , Ip % > 
+ = 0, (6) 
dr r 
would not be satisfied. On the other hand, if the material were able to 
flow within the radius r, in such a way that the thickness h were equal to 
hors/?, hg being the original thickness, the equation of equilibrium of a 
thin sheet of variable thickness, namely, 
d h 


(ho,) + —(o,—99) 0, (7) 
ar r 


would be satisfied. 
If the hole is expanded from a pin-hole to radius b and if the com- 
pressibility of the material be neglected, the volume of the metal contained 


between r = b and r = r, must be equal to argh). Thus 
Te 
. Cr hot: 
mrs ho —2-2(2xr dr), 
J fr 
b 
which yields b = 4r,. (8) 


Hence [A], _, = 2ho. (9) 
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Bethe’s model indicates a reason for the main feature of the observed | 


. : ;, is th 
plastic strains round the hole made in a sheet by a conical-headed punch aw 
or bullet, namely, the thickening of the plate close to the hole. When, | if 
however, actual plastic strain near the hole is considered in detail it will ee 
be seen that the arbitrary assumption that og = 0 and o, = —Y is | ‘itn 
inconsistent with any theory of plasticity in which the ratios of increments (11) 
of strain are related to corresponding ratios of stress differences. If pine 
op = o, = 0 and o, = —Y through the range b <r <1,, any theory of ae 
stress and strain in an isotropic plastic material would require that the x pa 
small increments in strain occurring while the radius of the hole was dista 
increased from 6b to b+-5b would be such that the radial strain e, is equal po 
to the normal strain e,. It is clear that this is not consistent with the actual iiead 


strain involved in the equation h = hyr,/r. In fact, for a small increment E 
66 in radius of the hole, Bethe’s model involves a radial displacement 5u, 

















cons 
where , stres 
du = (2—-—) db, pei 
( 5) dista 
, ee : had 
and the components of strain increment are 
tube 
Sep = 2b—r db lows 
r ob Ir 
dh = =4b b\ 5b ; a he 
$e, = — = —/{]1—-}— }. (10) 
ae r r} b The 
F 1 8b , 46°) 5b of ¢ 
i eo ee sym 
mm , ; , ; ; Ana 
Thus in Bethe’s analysis the ratios 5e,:5e9:5e, vary with radius while the y 
ratios of the principal stresses are constant. inf 
infix 
, _ : ad , Thu 
Problems in plasticity when the ratios of principal stress differ- 
. : eae stra 
ences are not constant at each point during the straining f 
8 fr 
process , 
m™. 3 eee. oe In t 
Che incompatibility, in any rational theory of plasticity, of Bethe’s ii 
hypothesis (og = o, = 0) with the strains represented by (10) must mean Ke 
that this hypothesis is incorrect. To improve the analysis it is necessary pha 
to use some hypothetical assumption or experimental result relating the the 
stress-difference ratios to the strain ratios. For small strains the simplest pid 
SS 
is that of Mises, namely, 
\ and 
O,— 99 __—«99 OZ sz D>, (11) also 
€,— 9 €g—e&, €,—e, , of t! 


This, though not quite an accurate representation of experimental results, 
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eved is the hypothesis most frequently used in discussing problems in plasticity. 
mech | tt will be used in the work which follows. 
hen, If the stress-difference ratios, namely, o,—og:0g—o,:0,—a,, are constant 
will at each elerhent of the material during the whole straining process, the 
r is correspondiig strain-difference ratios e,—é€g:eg—e,:e,—e, are constant and 
ents | (11) will apply if e,, eg, and e, are the total strains. If, however, the 
If stress-difference ratios change in the course of the deformation, (11) can 
y ot only apply if e,, eg, and e, are small increments of strain occurring during 
i the a portion of the deformation in which there is small change in the stress 
Was | distribution. Problems of plastic flow in which the stress ratios vary 
qual | during the deformation can therefore only be solved by following the 
tua strain-history of each element. 
—~ Eminent authorities have not always appreciated this point and have 
oo | consequently published erroneous solutions of problems in which the 
stresses and total plastic strains have been related as though the stress 
distribution had been constant during the deformation when in fact it 
| had not. The solutions of the problem of the straining of a thick-walled 
| tube beyond the elastic limit which are given by Nadait and by Soko- 
lowsky{ seem to me to be defective for this reason. 
In the following pages the solution is given to the problem presented by 
a hole in a sheet of plastic material which is expanded from a pin-hole. 
(29) The solution involves tracing the complete strain-history of each element 
of the sheet, but the analysis is much simplified by considerations of 
symmetry and similarity. 
4 Analysis of strain round an expanding radial hole in a sheet 
oa When a hole is enlarged the finite strain at any stage is made up of 
infinitesimal elements of strain which vary as the enlargement proceeds. 
For. Thus, when a small pin-hole in a plate is enlarged we must study the small 
ning strain produced in an element of the sheet which was originally at radius 
s from the pin-hole, when the hole enlarges from radius b to radius b+-8b. 
he’ In the more general case when the initial radius of the hole in the un- 
2 ; stretched sheet is not zero this is very difficult to analyse, but when the 
oa expansion starts from a small pin-hole it may be expected that the 
nerd configuration when the hole has radius 6, will be similar to that round 
ae the hole when its radius is b, except that the radii where any given thick- 
P hess occurs will be changed in the ratio b,/b,. Thus, if h is the thickness 
and w the radial displacement, it may be assumed that h/hy and u/b and 
(11) also the stresses are functions of s/b only, where hg is the initial thickness 
of the sheet, and s the initial radial distance of a particle. 
sults, {+ Plasticity (McGraw Hill, 1931), pp. 196-9. 


The Theory of Plasticity (Russian) (Moscow, 1946), Chapter ITI. 














108 G. I. TAYLOR 


To simplify matters I have assumed that the compressibility is so small 
that it may be neglected and the material taken as incompressible. The 
relationship between the small strain which occurs at any radius during 
the expansion of the hole through a small increase in radius from b to 
b+-86 can be understood by referring to Fig. 1. Here the ordinates repre- 
sent wu and the abscissae r. 


—tib b~ 























Fia. 1. 


The initial radial distances of the element, which at a subsequent stage 
in the opening-out of the hole is at radius 7, is related to u by the equation 


r= stu. (12) 


In Fig. 1, therefore, the displacement of a particle from its initial radius s 
is represented by a line drawn at 45° to the axes. In particular the 
displacement of the particles which were initially at the pin-point where 
the hole began is represented by the 45°-line OP,P,. The curved line 
P, AQ, represents the relationship between r and u which it is the object 
of the analysis to calculate. At a subsequent stage of the expansion, when 
the hole has expanded from radius b to radius b+-5b, the curve P, BCQ, 
representing displacement is similar to P, AQ ,, but with its linear dimen- 
sions increased in the ratio b+-5b:b; thus in Fig. 1, 


P,P) AC AD _ &b 


OP, AO ie 
so that AD = rdb/b. (13) 


If 5r is the change in r for a given particle of material when the hole 
expands from 6 to +46, dr is found by drawing the line AB at 45° to the 
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axes to meet the curve P, BCQ, in B. If 6b/b is small enough, the are CB 
may be taken as straight so that, if 7—a is the slope of CB to the axis, 
— = —tama. (14) 
or 
If 8 is the angle AOQ,, tan8 = u/r. From the geometry of the figure 
ABCD (Fig. 1) 
§r = AF BF = CEtana+DAtanfB = (DA—6r)tana+ DA tan B. 


| (15) 
Hence or [= n+ tan DA, (16) 
1-+-tan a 

‘—é6u/éer\ 8 
and, from (13), Sr (" titi. an (17) 
l—éu/er } 6b 


The radial strain component during the expansion of the hole from 6 to 
- € “¢r . . ‘a . ° o 
5b is — (Sr), and differentiating (17) with respect tor, keeping 5b constant, 
or - 
é | (l—u/r)r é?u] db 
( 


© (Sr) ms 
or" b 


. (18) 


1—du/dor)* or? 
Since the strain during expansion of the hole from 6 to b+-86 is propor- 
tional to 5b/b, it is convenient to define strain components e¢,, €g, and e, 
so that strains during the small enlargement 5b are «,5b/b, €95b/b, €,5b/b. 


With this definition 
- a2 
= [ptetel os 
(1—ou/or)* | or 
The tangential strain is simply 
b dr u/r—ou/er 


(20) 


€p _= 
6 ’ 7 ’ 
db r 1—ou/oer 


and the strain perpendicular to the sheet is 

€, —€,— €9. (21) 
The thickness h at any stage can be found simply from the equation of 
continuity: it is given by 


Bs (1—S)(1-3). (22) 
ho r or 


It can be verified that (21) is consistent with (22). 

These expressions for strain take simple forms when expressed in terms 

of a new independent variable é = r? and a new dependent variable 

7 = s* = (r—u)*. Making these transformations and writing 
dy : d*» 


=! ee 23 
P dé’ q dé?’ ( ) 
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(19), (20), and (21) become 


one 

ee Fie | 9 
€, + p Ep? (24) 

— oe 25 
€9 Ep’ (25) 

9 

ina: ana 9 

€z p ? (26) 


while (22) reduces to the simple form 
h/hy = p. (27) 
(26) can be deduced directly from (27). 
The stress-equilibrium equation for a thin sheet is 
o hia —0dp) 
—(ho,) + — - =— 0. (28) 
or r 
Two possible alternative forms for the strength condition might be 
considered: 
(a) Mohr’s stress criterion which may be written 


og—o, = Y if o@ is positive, i.e. tensile, 


(29) 
or —o,=Y if og is negative, i.e. compressive. 
(b) Mises’s condition which may be written, when o, = 0, 
o;-+0}—o90, = constant. (30) 


This reduces to —o, = constant if og = 0, and so is identical with Mohr’s 
in that case. 
If Bethe’s assumption that og = 0 combined with o, = constant is 
used, (28) leads to 
hr = constant = hgr), (31) 


where 7; is the outer boundary of the region of finite plastic strain. 


Substituting in (22) 
fs (3 _“\/1_ (32) 
r r ér)’ 


which gives on integration 
4(r—u)? = rr,+ constant. (33) 
Since u = 0 when r = 74, the constant is —4}(r,)? and 
u = r—,/{(2r—r9)r5}. (34) 
The inner boundary is where b = r = u, so that from (34) 
b = 4, (35) 


which is Bethe’s result if 7, is identified with his r,. 
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Introduction of the strain-ratio relationship 

Since the strain-ratio relationship (11) involves only differences between 
the principal stresses, we may assume without loss of generality that 
g,= 0. We shall assume further that the compressibility of the material 
i. small enough to be neglected so that «,4+¢9+¢«, = 0. In these circum- 


stances (11) may be written 


® . $4. (36) 





: , 
0, €, €. a€, 7 €9 


Substituting (27) and (36) in (28) the equilibrium condition reduces to 





d 0, [ €,—€ on 

2—(po,) +21 (6) — 0. (37) 
dé E \2e,+e9 

This equation must be used in conjunction with a strength criterion. 

Mohr’s criterion (a) will be used. In this case (37) assumes two different 

forms according as og is negative (i.e. compressive) or positive (i.e. tan- 


gential tension). These are: 


op negative, o, Y, so that (37) becomes 
5 " ) "= €9 ro . ‘ 
2q+ SS ) = 0; (38) 
S <¢, ! €9 
og positive, o,—ag Y, so that, from (36), 
Qe T €9 , 
o, = ———(—Y) 
€,— €9 
9 d 2€,+ €9 p ‘ 
and hence 2 —| pj ———} | += = 0. (39) 
dé €,—€9 é 





Substituting for «, and eg from (24) and (25) the resulting equations 


may be written: 


ag negative (tangential compression) 


>{ 4% ( 1) p n 
au— es ei Pe oe | 40 
a) 3)+e(—14 -" 
In this case, from (38), 
o_ _“*9tPis * 12 (41) 
é, 2q+p/é’ a, Pp’ 


and in terms of Mohr’s strength criterion the stresses are 


o, Y,: a= -Y (2), (42) 


O, 
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og positive (tangential tension) 


4n2q3 37? 3722 » In. xf 
{Ne} ee seo Bh Bag 


p* ép® Ep? ep] ' € Bp 
(43) 
3 
where w is written for i , 12. a 


The expressions for ¢g/e, and og/o, cannot be simplified by using the 
equation of equilibrium, and the full expressions derived from (24), (25), 
and (36) must be used, namely, 

- 1—n/Ep Se ___ ng an | (44) 
« —1+2nq/p*+n/fp 9, = 4nqt+npl——p> 


and in terms of Mohr’s condition o,—o, = —Y the stresses are now 


= —Tit—oio n= —¥)/(-2). 
0, 0, 


It will be seen that (43) is an ordinary differential equation of the third 
order and first degree while (40) is of the second order and second degree. 
The reason for this difference lies in the form of Mohr’s strength condition. 
When ay is positive three boundary conditions can be assigned at any given 
value of € (i.e. of r). These might, for instance, be w/r, h/hy, and o, which 
can be transformed directly in assigned values of g, p, and y. When o, 
is negative o, cannot be assigned arbitrarily; it is in fact constant. Thus 
only p and 7 can be assigned arbitrarily. 


Boundary condition at the elastic-plastic boundary 


The elastic stresses due to radial displacement in an infinite sheet are 





—9, = o9 = $Yr4/?", (46) 
where r, is the radius at which o,—og = —Y. The corresponding small 
radial displacement is 

l ar -2 
“u= —_ = (=)2. (47) 


where # is Young’s modulus and m is Poisson’s ratio. In the present 
investigation compressibility will be neglected and we will take m = }. 
In the elastic region, therefore, where u/r is small compared with unity, 
u\2 3 VY 72 
y= r{1—-) =¢ _3¥n\ (48) 
r ’ 2H €é 
At the inner boundary of the elastic region therefore 


we 


dn 
SS as GS i. == == 
Pp dé ; q = 0, i] e( 








At tl 
crite! 
throu 
in th 
hY. | 
boun 
be co 
09 be 
the p 
Stra 

In 
foun 


wher 


and 


Whe 


Subs 


pare 


wher 


Subs 


Neg 
do x 


The 
b/(u 


C fo 


(44) 


hird 
Tee, 
tion. 
iven 
hich 
n 0p 
Thus 


are 
(46) 


mall 








A CIRCULAR HOLE IN A THIN PLASTIC SHEET 113 


At the outer boundary of the plastic region, since o, is positive, Mohr’s 
criterion ensures that o,—og = —Y’. Since o, is necessarily continuous 
through 7 = r,, and it is assumed that o,—og = —Y at the elastic limit 
in the elastic region, og must be continuous through r = r, and equal to 
iY. It is important to notice the reason why og is continuous at the plastic 
boundary in this case, because it is not necessary in general that og shall 
be continuous when Mohr’s criterion is used. It will be shown in fact that 
og becomes discontinuous on the circle r = r, within the plastic region at 
the point where og = 0. 


Strains and displacements when 7, >r >r, 
In the region within the circle r = r,, where og is positive, it will be 
found that the strains are small, being of order Y/H. Assuming that 


n = &(l—an,) and p= 1+ap,, 
where x = 3Y/2E, 
dy d d - 
p= = = 1—an,—a€ 3 so that p, = —n—8 ae (50) 
Ip dp d d? 
a ; aj _ _, —— g@N dl 5 | ; x 
- a a” Ee ae +E aes _ 
When og is positive 
a ee ee 
o, y(= 7 - -w(= rnlé+4ng r), (52) 
 &—€ —p+nlé+ngq/p 


Substituting from (50) and (51) in (52) and neglecting terms in a? com- 
pared with those containing «, 


a= (e422) 


| if = 
where ab z and yf’ = a 
Substituting this in (39), 
d pus P a 
an as 54 
dé (7 a " 3é os 


Neglecting terms which contain «a as a factor compared with those that 
do not, » may be taken as 1-0 and (54) may then be integrated giving 


us 


—+1]né = constant. (55) 
b+ Eus 7 7 
The boundary condition at r = r, is ¢, = —og = —}Y so that from (53) 
W (bt Ep’) 1. The constant in (55) is therefore —1+4lnrj. Writing 
¢ for In(é rz), (565) becomes 
dis 6+ 
+62) = ©, (56) 
d¢ 3-+¢ 


5092.1 


I 
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The integral of (56) is 


In’+lC+3In(3+2) = constant. (57) 
Throughout the whole of the elastic region —o, = og and o, = 0, so that, 


at the elastic boundary, p = 1; hence, from (50), 








| Ko 
m+éh| = 0, (58) 
| dé on 
and from the definition of ,, (49) shows that [7]. = 1. Hence from 
(58), [Ee]. » = —1; (57) may therefore be written 
#é(1+4Inéry?)3 = —1, (59) 
. dy, «i : 
and since 4 = dé and € = éry;* this becomes 
dn = = — (60) 
dl (1+3¢)° 
(60) may be integrated, giving 
3 1 ; 
ne 20440)" 2° 8h) 
The equation for the thickness is 
h dy dy 1 1+2 : 
— = = f— ——) = l-+al—-— ~om (62) 
hy dé ant dé “12 204403 
The displacement is 
BS = £i— y} = danr = bar[ $(1—4f)-?—}]. (63) 
' ae 6+ ld 
The stress can be found by substituting from (56) — - > for — we i.e. for 
? 3d+¢ yb d 
a , in (53). It is found that 
a 
o, = —4Y(1—f) = —4Y(1—2In(r/r,)}, (64) 
and hence og = Y—o, = ay(1 +2In ‘), (65) 
r, 


This is the well-known result which can be obtained without considering 
the strains and displacements, if it is assumed that the thickness of the 
plate does not vary. 


Values when o, = 0 

The radius r at which og = 0 is from (50) rz = r,e~! = 0-606r, and 
corresponds to ¢ = —1. 

Though the stress distribution in the range 7, > 7 > 7, is identical with 
that found in Bethe’s investigation, the displacements and strains are not 








the 
giv 


Put 
hy} 


The 


Boi 


A 
the 


Sub 
resv 


Nei 
the 


of o 


a pr 
gets 
hav 








Om 


60 


(63 


. for 


(64 
(65 


Ting 
’ the 


and 


with 


> not 








A CIRCULAR HOLE IN A THIN PLASTIC SHEET 115 


the same. For the case when Poisson’s ratio is } Bethe’s method would 
give for the displacement when og = 0 


9 
ry Ler d > Fal 
Us tx = far,(2-718). (66) 
2 tr 
Putting ¢ | in (63) the displacement according to the present strain 
hypothesis is 
xr. 3 l 23 wien —_ 
Us 2| — = tar,|—] = 4ar,(2°875). (67) 
F 2 |2(2)? 2) (38 7 


The displacement is in fact about 6 per cent. greater than that calculated 
on Bethe’s strain hypothesis. 

Putting ¢ 1 in (62) the value of h/hy at r = r, is 1+ 4a, and from 
(51) and (56) 


i dil\ q 6+¢\ __ ¢ \. 
q (29 o- zi) :— (2-5) ss —ob(s 


hence, from (57), ge = x6 
so that when ¢ —l, 
qg& = —j6e- (68) 


Boundary values at op = 0 


At the circle r = r,, where og = 0, p and 7 are continuous. Just inside 


the circle therefore, where og is negative, 
1 o« h , " 
: Sy and p= — = 1+4a. (69) 
hy 


l . . 
x1) l 
S 
When og is negative q is determined by (40) when »/é and p are given. 
0 1 A 7 ; 
Substituting in (40) from (69), the values of gé found by solving the 
resulting quadratic equation are (neglecting terms in a?) 
gf 1 Sa and g&é= +%a. (70) 
Neither of these values is the same as q = —27a/16, the value just outside 
the boundary, so that q is not continuous at r = rp. 
The two possible values of g just inside the radius r = r, are therefore 
of opposite signs. Since 
dp 1 dh ] dh 
q a —_- = 7 BP 
dé hyd& = 2hgr dr’ 
a positive value of g would correspond to a condition in which the plate 
gets thinner towards the centre, a configuration which does not appear to 
have physical significance. The negative value of q for which gé = —}—%3a 
will therefore be taken as the boundary value to be used in continuing the 
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solution inward from r = r. It will be noticed that ¢,, and consequently 
og, are discontinuous as well as g. This discontinuity arises from the form 
of Mohr’s criterion. It would not occur if von Mises’s criterion had been 
used. 


Discontinuity in «, and eg 


Substituting 
a= —t-He, palit, = 1-H 
in (24) and (25), it is found that 
=a, =i Ts, 
and, substituting these in (36), 
G, ~ 
2 = 4—fa 
Cd. 


When « is small, i.e. when £/Y is small, we may neglect « and take as the 
boundary condition at 7 = r, for calculating the stresses and displacements 
when r < 1’, the values 


p=1, wf=1, €=-} (71) 
and the stresses are o, = —Y, og = —}4Y. 


Thus the stress og suddenly changes from 0 to a compressive stress of 
—+Y at the radius r = 7. 


Calculation of stress and strain when r < 7, 

To calculate the distribution of stress and plastic strain inside the 
radius r= 73, (40) must be solved step by step using the boundary 
values (71). 

If 5€ is the magnitude of a small step, the corresponding changes p and 7 
may be taken as 

dn = pdé+43q(dE)? ) (72) 

bp = qdé J " 
After calculating the values of p and 7 at the end of each step these values 
are inserted in (40) and the resulting quadratic for q is solved, the root 
which derives by continuous variation of 7, and p from qr? = —}, being 
chosen in each case. 

The results of applying this process are given in Table 1 and are shown 
graphically in Fig. 2. Values of the principal variables £/r3 and »/r} are 
given in cols. 1 and 2, Table 1. Values of p and —qr3 are given in cols. 3 
and 4. Using these values of p, g, and €, values of ag/o, calculated from 
(41) are given in col. 8, and the corresponding values of o,/Y and o/Y in 
cols. 9 and 10. It will be seen that op, which begins as a compressive 
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stress equal to half the radial stress at the outer limit of the region of 
finite plastic flow, rapidly decreases till when ¢/r3 = 0-35 it becomes zero, 
and if the process is carried farther, using (41), og becomes a tension. 











r 
b 
Fia. 2. 


When ¢/r? = 0-30, for instance, the calculated value of o/c, is —0-124. 
For values of £/r3 less than 0-35, therefore, the alternative form (43) of the 
equilibrium equation must be used. 














TABLE 1 
I 2 4 z 6 7 5 } 9 10 
ely 7 p yr2 ord r/b u/b 09/c, o,/Y | o/¥ 
; I I 5 2°21 fc) 0°50 are) | —o'50 |) 
090 0899 1°025 305 27096 | ‘oor | | —1°o —0°47 
080 79 1°OS55 o°381 1'°978 | +008 +0°440 | —I —0°44 
075 741 1°075 0°431 1-915 "012 | +0°397 | —I°0 —0o°40 |j Us 
o"7 ” 1°096 0°493 1°850 “O19 0°370 I°o —0°37 $ > 
0°65 0°632 I°r21 0°566 1°782 | *o2 0°343 I°o —O34 ll ¢ 3 
; STE I*149 0°66 1-712 | 035 +o310|—-ro | —o31 |f 3 ® 
050 457 1°257 2-998 | 1°563 | 7070 «| +0212} —1r0 =| —ov2I Ee 
on45 392 1°317 1°24 1°483 | ‘100 | +0152 | —1-0 —o15 || & ~ 
0°40 325 1°379 1°583 1°398 | +138 | +0°082 | —1-o -0'08 || 
0°35 254 I°450 2070 1-307 | +194 | +0°000 | —1-0 ° 
0 179 1°554 2910 1-210 *276 =| —o-124 ro) =| +012 |) 
35 0254 1°450 2°O70 26°5 1308 | +194 a 
0°30 178 1°587 3°397 57°3 I*210 278 0092 | —0°916 | +0°084| | 3 3 
7 129 I°715 5117 108-0 1°149 *356 0168 0°857 | +0°143|| ~ S 
1 7 ‘917 8357 81-0 1-083 -478 0°328 | —0°753 | +0°247/ > =. z 
34 2°140 | 13°98 )34°0 1037 "630 =| —0°569 0°638 | +0°362| : = 
326 | 23°22 51900 1013 771 0°739 | —0°576 | +0-424/| & © 
205 2°61 1°00 =| ‘I-00 -1*000 | —0°500 | +0500) | ~ 





Since o, is continuous and equal to —Y at ¢/r} = 0-35 and og = 0 when 
é/r? is just greater than 0-35, while Mohr’s criterion ensures that 
o,— 09 = —Y 
when og is positive, it seems that og = 0 when é/r2 is just less than 0-35. 
Since both o, and og are therefore in this case continuous through the 
tadius where og changes sign, ¢, and eg are also continuous. Hence from 
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(24) g is continuous. The values of 7, p, and q at €/r} = 0-35 can therefore 
be inserted in (43) and the value of w = d°7/d& at €/r3 = 0-35 determined. 

The changes in 7, p, and q during the first step 5€ in the new region are 
calculated using the formulae 


by = pdt+3q(df)?+ gw(dé)?, 

dp = qd€+ 4u(dEé), (73) 

dq = woé. 
Values of »/r2, p, —qr3, and wr} found in this way are given in the lower 
part of Table 1, corresponding to 0-35 > &/r§ > 0-205. Values of o9/c,, 
o,/Y, and o9/Y from (44) and (45) are given in cols. 8, 9, 10 of Table 1. 


Conditions at edge of hole 

It will be seen in Table 1 that as €/r} decreases to 0-21, —qr3 and wr} 
are rising very rapidly. A study of the values of the terms in (43) reveals 
that by the time €/r} = 0-21 is reached, one term on the right-hand side 
of the equation and one on the left-hand side are larger than any other 
terms. The limiting form of the equation when 7 is small is in fact 


3w 
i — —29. (74) 
p 
This equation can be integrated twice, thus 
— — An ~% | 
(75) 
and p? = B—6Ax} J 
A and B being the two constants of integration. To determine the 
constants the values of /r} = 0-012, p = 2-326, gr? = —23-22 can be 
used at € = 0-21. The resulting values of A and B are 
A — 1-217rz8, B = 7:07. (76) 


The limiting value of p when 7 = 0 is therefore 
p = V7-07 = 2°66. (77) 


This is the limiting value of h/h, at the edge of the hole and may be 
compared with Bethe’s value 2-0. It is not very different from the value 
at €/r{ = 0-21; to find the limiting value of € therefore it is sufficient to 
take p as constant and equal to 2-61 in the interval during which 7 
decreases from 0-012 to 0. Thus the limiting value of €/r2 corresponding 
to the edge of the hole is 


0-012 78 
é. 6 —* — 0-21 — = 0-205. (78) 


2-66 
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re i ._ (radius of finite plastic deformation\ . 
[The ratio : : is : 
dd. radius of hole 
ire 
Ys ] = 
- = ce 2°21 (79) 
0 V0°205 
This may be compared with Bethe’s value 2-0. 
3) Substituting the approximate limiting forms of p and q from (75) in 


44), the limiting form for a/c, is 


r 


- ; 7-30ntr; § 
- lim (7°) lim —— 1 "2  ¥ (80) 


= i P —? 
C.. 7-0 0, n—0 1°O4 2-43y'r5 





This tends to the value —1 as indicated in the last figure of col. 8, and 
the corresponding values of o, and og are therefore —0-5Y and +0-5Y. 
It will be noticed that the stress at the internal boundary could have 


prt | been predicted a priori if it had been possible to assume that h/hg is finite 
ale at r= b, because clearly the total amounts of strain in the tangential 
de and radial directions are both infinite at a hole which has been enlarged 
er from a pin-hole. Thus the state of strain at the hole is such that symmetry 
alone must ensure that o, is exactly half-way between o, and og. Since 
, = 0, ¢, og. Similar considerations can be used to understand why 
/4) the stress at points just inside the boundary r = r, corresponds with 
| (a9/0,) _(-5, for at the edge of the region of finite plastic displacement, 
| where the radial displacement is zero, eg = 0. Thus «, = —e, and og must 
| therefore be exactly half-way between o, and o,. Hence, since o, = 0, 
75) 19 = t0,. 
he Expressions in terms of radius of hole 
be | The radial variable is expressed in terms of the radius of the plastic 
region. To express the results in terms of 6, it is necessary to tabulate 
+6) b = 2-21,é/r,. These values are given in col. 6, Table 1. The displace- 
ments u/b = 2-21(.€—~yn)r,' are tabulated in col. 7. 
The radial displacement is shown graphically in Fig. 2 which may be 
77) ompared with the diagrammatic sketch, Fig. 1. 
be Comparison with Bethe’s results 
lue The thickness ratio h/hy = p is shown in Fig. 3 and Bethe’s values, 
to namely, h/h, = 2b/r, are also shown. It will be seen that the main differ- 
L 1 ences are that the present calculation shows the ‘crater’ extending farther 
ing radially than Bethe’s and at the same time the ‘crater’ is much steeper 
‘lose to the hole. From (75) it will be seen that q> —® as » > 0 and, 
78) since q 1 dh the equati s indicate the > edg j > ‘crater’ is 
, juations indicate that the edge of the ‘crater’ is 


2h, r dr 
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a thin knife-edge. This deduction, however, is based on the assumption 





that the strains and stresses are uniform through the thickness of the plate. | 
It is evident that this assumption ceases to be a good approximation when 


dh . F ; ’ 
— > —oo. At first sight it might be thought that the extra thickness at 


dr 
= 6 above Bethe’s 2h means that the work done in expanding the hole 
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is greater according to the present calculations than in Bethe’s calcula- 
tions, but this is the reverse of the truth, for the radial stress at the hole 
is only —}Y instead of Bethe’s —Y. In fact the work done in expanding | 
to a given radius is 

mb*ho(4Y)(2-66) = 1-337b%h, Y (81) 


while if Bethe’s strain assumption is used it is 27b*h, Y. 

Fig. 4 shows the distribution of stress. This is of course very different 
from Bethe’s, the most striking difference being that the present calcula- | 
tions predict a state of tangential tension in a ring which extends to 
30 per cent. of the radius of the hole from its edge and a tangential 
compression from that point to the edge of the region of large plastic 
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distortion. In the plastic region r,; >r >r,, where small strains com- 
parable with the elastic strains occur, the stress is as calculated by Bethe, 
ie. there is a tangential tension. In this connexion it may be noticed that 
in comparing calculations of this kind with the behaviour of real materials, 
a metal which experiences considerable hardening with cold work might 
give results differing widely from the above theory. 


Experimental work 


An attempt was made to produce the distribution of strain contemplated 
in the foregoing analysis. A series of tapering sections of a cone of semi- 
vertical angle 14° were made in mild steel. The smallest of these, which 
terminated in the vertex, was mounted in the chuck of a lathe and a flat 
sheet of lead was mounted on a slide rest so that its plane was perpendicular 
to the lathe bed. A small hole was bored in the lead sheet and the slide 
rest was moved parallel to the axis of rotation of the cone so that the 
vertex entered the hole and exerted radial pressure on its circumference. 
The cone was rotated while the hole was being enlarged so that the 
direction of the friction should be in the tangential direction rather than 
perpendicular to the plane of the lead sheet. ‘Oildag’ was used as a high- 
pressure lubricant between the sheet and the cone. 

On reaching the base of the smallest section the lead sheet was withdrawn 
and the next largest section of the cone was fitted in the chuck. The hole 
was then enlarged to the largest diameter of this section. This process was 
repeated with successive sections, or broaching tools, till the hole was 
2-45 cm. in diameter. 

It was found that the hole could be enlarged symmetrically without 
bending the plate till its diameter was between 7 and 10 times the thickness 
ifthe sheet, but that at about this stage the sheet always bent out of its 
plane. This was not due to a sideways pressure, for the sheet sometimes 
vent towards the thick end of the broaching tool and sometimes away 
irom it. It appears that the sheet is unstable and that an alternative form 

if deformation occurs in which the ‘crater’ develops into a short length 
of tube on one side or other of the sheet. 

This tube is joined to the flat sheet by a region where the sheet is bent 
into the form of a curved fillet. The wall of the tube varies in thickness 
irom that of the sheet at the fillet to zero at the outer edge. This edge 
contains the particles which were originally at the point where the vertex 
of the cone entered the sheet. 

Sections of deformed sheets are shown in Fig. 5. These photographs 
were obtained by sawing the sheets in two in planes passing through their 
ixes of symmetry and grinding down the surfaces so obtained till the true 
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sections were revealed. The upper photograph shows the hole when its 
radius is 2-7 times the thickness of the sheet, while in the lower photograph 
the hole is about 10 times the thickness. It is clear that in the unsymmetri- 
cal case (lower photograph, Fig. 5) the material of the tube-shaped crater 
is expanding under the influence of a tangential tension Y and that the 
stress perpendicular to the plane of the sheet is zero. The radial stress 
varies from zero at the outer surface to approximately Yt/b at its inner 
surface. Here 6 is the radius of the hole and ¢ the thickness of the wall 
of the tubular crater. Since ¢ is less than hy the radial stress is negligible 
when b/h, is large. The stress is therefore approximately uni-axial. The 
contraction of the material in the radial direction should therefore be 
equal to the contraction in the direction perpendicular to the sheet. This 
type of deformation would give rise to a distribution of thickness ¢ in the 
wall of the tubular crater for which t = ho,/(y/D), where y is the distance 
from the plane of the lip of the crater and D is the height of this plane 
above the plane of the sheet. The initial radius, s, of the ring of particles 
which are at the point in the crater defined by y is given by 


Ss y 89 
b [5 (82) 


Since the material is assumed to be incompressible 


D 
1b*h, = 2mb } Y dy, (83) 
0 
so that D = 3b (84) 
and t h y _- = 1-15h, ly (85) 
"A 0-756 °/ b 


The main features of the approximate theory represented by (84) and 
85) are found, at any rate qualitatively, in the experiment represented 
in the lower photograph of Fig. 5. Since this unsymmetrical mode of 
deformation and also the symmetrical deformation are both theoretically 
possible it is of interest to compare the work done in opening out the 
hole to a given radius b by the two modes. 

The work done in expanding the ring of material which was originally 
contained between radii s and s+5s is 


: is ij 
2a} hos d5slin : 
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so that the total work done is 


b 
w = 2h, Y | sin2ds = drb*hy Y. (86) 
8 
0 


Comparing (86) with (81), it will be seen that the symmetrical mode 
requires 2-6 times as much work as the unsymmetrical mode for the same 
radius of hole. This fact seems to explain why it is the unsymmetrical 
type of deformation which actually occurs. 
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